65 MyLength_(map.NumMyPoints()),
66 GlobalLength_(map.NumGlobalPoints64()),
67 NumVectors_(numVectors),
68 UserAllocated_(false),
69 ConstantStride_(true),
70 Stride_(map.NumMyPoints()),
90 MyLength_(Source.MyLength_),
91 GlobalLength_(Source.GlobalLength_),
92 NumVectors_(Source.NumVectors_),
93 UserAllocated_(false),
94 ConstantStride_(true),
101 double ** Source_Pointers = Source.
Pointers();
112 double *
A,
int MyLDA,
int numVectors)
117 MyLength_(map.NumMyPoints()),
118 GlobalLength_(map.NumGlobalPoints64()),
119 NumVectors_(numVectors),
120 UserAllocated_(false),
121 ConstantStride_(true),
122 Stride_(map.NumMyPoints()),
142 double **ArrayOfPointers,
int numVectors)
147 MyLength_(map.NumMyPoints()),
148 GlobalLength_(map.NumGlobalPoints64()),
149 NumVectors_(numVectors),
150 UserAllocated_(false),
151 ConstantStride_(true),
152 Stride_(map.NumMyPoints()),
173 int *Indices,
int numVectors)
178 MyLength_(Source.MyLength_),
179 GlobalLength_(Source.GlobalLength_),
180 NumVectors_(numVectors),
181 UserAllocated_(false),
182 ConstantStride_(true),
191 double ** Source_Pointers = Source.
Pointers();
204 int StartIndex,
int numVectors)
209 MyLength_(Source.MyLength_),
210 GlobalLength_(Source.GlobalLength_),
211 NumVectors_(numVectors),
212 UserAllocated_(false),
213 ConstantStride_(true),
222 double ** Source_Pointers = Source.
Pointers();
263 int randval = rand();
267 int locrandval = randval;
295 #ifdef EPETRA_HAVE_OMP 296 #pragma omp parallel for default(none) shared(to,from) 297 for (
int j=0; j<myLength; j++) to[j] = from[j];
299 memcpy(to, from, myLength*
sizeof(
double));
317 int randval = rand();
321 int locrandval = randval;
362 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 366 EPETRA_CHK_ERR(ChangeGlobalValue<int>(GlobalRow, 0, VectorIndex, ScalarValue,
false));
371 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 375 EPETRA_CHK_ERR(ChangeGlobalValue<long long>(GlobalRow, 0, VectorIndex, ScalarValue,
false));
380 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 382 int VectorIndex,
double ScalarValue) {
384 EPETRA_CHK_ERR(ChangeGlobalValue<int>(GlobalBlockRow, BlockRowOffset, VectorIndex, ScalarValue,
false));
389 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 391 int VectorIndex,
double ScalarValue) {
393 EPETRA_CHK_ERR(ChangeGlobalValue<long long>(GlobalBlockRow, BlockRowOffset, VectorIndex, ScalarValue,
false));
398 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 402 EPETRA_CHK_ERR(ChangeGlobalValue<int>(GlobalRow, 0, VectorIndex, ScalarValue,
true));
407 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 411 EPETRA_CHK_ERR(ChangeGlobalValue<long long>(GlobalRow, 0, VectorIndex, ScalarValue,
true));
416 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 418 int VectorIndex,
double ScalarValue) {
420 EPETRA_CHK_ERR(ChangeGlobalValue<int>(GlobalBlockRow, BlockRowOffset, VectorIndex, ScalarValue,
true));
425 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 427 int VectorIndex,
double ScalarValue) {
429 EPETRA_CHK_ERR(ChangeGlobalValue<long long>(GlobalBlockRow, BlockRowOffset, VectorIndex, ScalarValue,
true));
442 int VectorIndex,
double ScalarValue) {
455 int VectorIndex,
double ScalarValue) {
461 template<
typename int_type>
463 int VectorIndex,
double ScalarValue,
bool SumInto) {
465 if(!
Map().
template GlobalIndicesIsType<int_type>())
466 throw ReportError(
"Epetra_MultiVector::ChangeGlobalValues mismatch between argument types (int/long long) and map type.", -1);
474 int VectorIndex,
double ScalarValue,
bool SumInto) {
478 if (BlockRowOffset<0 || BlockRowOffset>=
Map().ElementSize(MyBlockRow))
EPETRA_CHK_ERR(-2);
483 Pointers_[VectorIndex][entry+BlockRowOffset] += ScalarValue;
485 Pointers_[VectorIndex][entry+BlockRowOffset] = ScalarValue;
507 #ifdef EPETRA_HAVE_OMP_NONASSOCIATIVE 508 #pragma omp parallel for default(none) 510 for(
int j = 0; j < myLength; j++)
528 double * to =
A + i*MyLDA;
529 for (
int j=0; j<myLength; j++) *to++ = *from++;
561 double * to = ArrayOfPointers[i];
562 memcpy(to, from, myLength*
sizeof(
double));
602 #ifdef EPETRA_HAVE_OMP 603 #pragma omp parallel for default(none) shared(ScalarConstant) 605 for (
int j=0; j<myLength; j++) to[j] = ScalarConstant;
623 int *PermuteFromLIDs,
631 double **From =
A.Pointers();
635 int * ToFirstPointInElementList = 0;
636 int * FromFirstPointInElementList = 0;
637 int * FromElementSizeList = 0;
641 if (!ConstantElementSize) {
643 FromFirstPointInElementList =
A.Map().FirstPointInElementList();
644 FromElementSizeList =
A.Map().ElementSizeList();
654 if (MaxElementSize==1) {
656 NumSameEntries = NumSameIDs;
658 else if (ConstantElementSize) {
660 NumSameEntries = NumSameIDs * MaxElementSize;
664 NumSameEntries = FromFirstPointInElementList[NumSameIDs];
668 if (To==From) NumSameEntries = 0;
673 for (
int i=0; i < numVectors; i++) {
674 if (CombineMode==
Epetra_AddLocalAlso)
for (
int j=0; j<NumSameEntries; j++) To[i][j] += From[i][j];
675 else for (
int j=0; j<NumSameEntries; j++) To[i][j] = From[i][j];
679 if (NumPermuteIDs>0) {
685 if (CombineMode==
Epetra_AddLocalAlso)
for (
int j=0; j<NumPermuteIDs; j++) To[0][PermuteToLIDs[j]] += From[0][PermuteFromLIDs[j]];
686 else for (
int j=0; j<NumPermuteIDs; j++) To[0][PermuteToLIDs[j]] = From[0][PermuteFromLIDs[j]];
689 for (
int j=0; j<NumPermuteIDs; j++) {
690 jj = PermuteToLIDs[j];
691 jjj = PermuteFromLIDs[j];
692 if (CombineMode==
Epetra_AddLocalAlso)
for (
int i=0; i<numVectors; i++) To[i][jj] += From[i][jjj];
693 else for (
int i=0; i<numVectors; i++) To[i][jj] = From[i][jjj];
700 for (
int j=0; j<NumPermuteIDs; j++) {
701 jj = MaxElementSize*PermuteToLIDs[j];
702 jjj = MaxElementSize*PermuteFromLIDs[j];
703 if (CombineMode==
Epetra_AddLocalAlso)
for (
int i=0; i<numVectors; i++)
for (k=0; k<MaxElementSize; k++) To[i][jj+k] += From[i][jjj+k];
704 else for(
int i=0; i<numVectors; i++)
for (k=0; k<MaxElementSize; k++) To[i][jj+k] = From[i][jjj+k];
711 for (
int j=0; j<NumPermuteIDs; j++) {
712 jj = ToFirstPointInElementList[PermuteToLIDs[j]];
713 jjj = FromFirstPointInElementList[PermuteFromLIDs[j]];
714 int ElementSize = FromElementSizeList[PermuteFromLIDs[j]];
715 if (CombineMode==
Epetra_AddLocalAlso)
for (
int i=0; i<numVectors; i++)
for (k=0; k<ElementSize; k++) To[i][jj+k] += From[i][jjj+k];
716 else for (
int i=0; i<numVectors; i++)
for (k=0; k<ElementSize; k++) To[i][jj+k] = From[i][jjj+k];
741 double **From =
A.Pointers();
746 int * FromFirstPointInElementList = 0;
747 int * FromElementSizeList = 0;
749 if (!ConstantElementSize) {
750 FromFirstPointInElementList =
A.Map().FirstPointInElementList();
751 FromElementSizeList =
A.Map().ElementSizeList();
754 double * DoubleExports = 0;
756 SizeOfPacket = numVectors*MaxElementSize*(int)
sizeof(
double);
758 if(SizeOfPacket*NumExportIDs>LenExports) {
759 if (LenExports>0)
delete [] Exports;
760 LenExports = SizeOfPacket*NumExportIDs;
761 DoubleExports =
new double[numVectors*MaxElementSize*NumExportIDs];
762 Exports = (
char *) DoubleExports;
767 if (NumExportIDs>0) {
768 ptr = (
double *) Exports;
771 if (MaxElementSize==1) {
774 for (
int j=0; j<NumExportIDs; j++)
775 *ptr++ = From[0][ExportLIDs[j]];
778 for (
int j=0; j<NumExportIDs; j++) {
780 for (
int i=0; i<numVectors; i++)
781 *ptr++ = From[i][jj];
787 else if (ConstantElementSize) {
789 for (
int j=0; j<NumExportIDs; j++) {
790 jj = MaxElementSize*ExportLIDs[j];
791 for (
int i=0; i<numVectors; i++)
792 for (k=0; k<MaxElementSize; k++)
793 *ptr++ = From[i][jj+k];
800 int thisSizeOfPacket = numVectors*MaxElementSize;
801 for (
int j=0; j<NumExportIDs; j++) {
802 ptr = (
double *) Exports + j*thisSizeOfPacket;
803 jj = FromFirstPointInElementList[ExportLIDs[j]];
804 int ElementSize = FromElementSizeList[ExportLIDs[j]];
805 for (
int i=0; i<numVectors; i++)
806 for (k=0; k<ElementSize; k++)
807 *ptr++ = From[i][jj+k];
833 if( CombineMode !=
Add 834 && CombineMode !=
Zero 841 && CombineMode !=
AbsMax )
844 if (NumImportIDs<=0)
return(0);
851 int * ToFirstPointInElementList = 0;
852 int * ToElementSizeList = 0;
854 if (!ConstantElementSize) {
862 ptr = (
double *) Imports;
865 if (MaxElementSize==1) {
868 if (CombineMode==
InsertAdd)
for (
int j=0; j<NumImportIDs; j++) To[0][ImportLIDs[j]] = 0.0;
869 if (CombineMode==
Add || CombineMode==
InsertAdd)
for (
int j=0; j<NumImportIDs; j++) To[0][ImportLIDs[j]] += *ptr++;
870 else if(CombineMode==
Insert)
for (
int j=0; j<NumImportIDs; j++) To[0][ImportLIDs[j]] = *ptr++;
871 else if(CombineMode==
AbsMax)
for (
int j=0; j<NumImportIDs; j++) {To[0][ImportLIDs[j]] =
EPETRA_MAX( To[0][ImportLIDs[j]],std::abs(*ptr)); ptr++;}
872 else if(CombineMode==
AbsMin)
for (
int j=0; j<NumImportIDs; j++) {To[0][ImportLIDs[j]] =
EPETRA_MIN( To[0][ImportLIDs[j]],std::abs(*ptr)); ptr++;}
873 else if(CombineMode==
Epetra_Max)
for (
int j=0; j<NumImportIDs; j++) {To[0][ImportLIDs[j]] =
EPETRA_MAX( To[0][ImportLIDs[j]],*ptr); ptr++;}
874 else if(CombineMode==
Epetra_Min)
for (
int j=0; j<NumImportIDs; j++) {To[0][ImportLIDs[j]] =
EPETRA_MIN( To[0][ImportLIDs[j]],*ptr); ptr++;}
875 else if(CombineMode==
Average)
for (
int j=0; j<NumImportIDs; j++) {To[0][ImportLIDs[j]] += *ptr++; To[0][ImportLIDs[j]] *= 0.5;}
901 for (
int j=0; j<NumImportIDs; j++) {
903 for (
int i=0; i<numVectors; i++) {
904 if (CombineMode==
InsertAdd) To[i][jj] = 0.0;
905 if (CombineMode==
Add || CombineMode==
InsertAdd) To[i][jj] += *ptr++;
906 else if (CombineMode==
Insert) To[i][jj] = *ptr++;
907 else if (CombineMode==
AbsMax) {To[i][jj] =
EPETRA_MAX( To[i][jj], std::abs(*ptr)); ptr++; }
908 else if (CombineMode==
AbsMin) {To[i][jj] =
EPETRA_MIN( To[i][jj], std::abs(*ptr)); ptr++; }
911 else if (CombineMode==
Average){To[i][jj] += *ptr++; To[i][jj] *= 0.5;}}
965 else if (ConstantElementSize) {
967 for (
int j=0; j<NumImportIDs; j++) {
968 jj = MaxElementSize*ImportLIDs[j];
969 for (
int i=0; i<numVectors; i++) {
970 if (CombineMode==
InsertAdd)
for (k=0; k<MaxElementSize; k++) To[i][jj+k] = 0.0;
971 if (CombineMode==
Add || CombineMode==
InsertAdd)
for (k=0; k<MaxElementSize; k++) To[i][jj+k] += *ptr++;
972 else if (CombineMode==
Insert)
for (k=0; k<MaxElementSize; k++) To[i][jj+k] = *ptr++;
973 else if (CombineMode==
AbsMax) {
for (k=0; k<MaxElementSize; k++) { To[i][jj+k] =
EPETRA_MAX( To[i][jj+k], std::abs(*ptr)); ptr++; }}
974 else if (CombineMode==
AbsMin) {
for (k=0; k<MaxElementSize; k++) { To[i][jj+k] =
EPETRA_MIN( To[i][jj+k], std::abs(*ptr)); ptr++; }}
975 else if (CombineMode==
Epetra_Max) {
for (k=0; k<MaxElementSize; k++) { To[i][jj+k] =
EPETRA_MAX( To[i][jj+k], *ptr); ptr++; }}
976 else if (CombineMode==
Epetra_Min) {
for (k=0; k<MaxElementSize; k++) { To[i][jj+k] =
EPETRA_MIN( To[i][jj+k], *ptr); ptr++; }}
977 else if (CombineMode==
Average) {
for (k=0; k<MaxElementSize; k++) { To[i][jj+k] += *ptr++; To[i][jj+k] *= 0.5;}}
1037 int thisSizeOfPacket = numVectors*MaxElementSize;
1039 for (
int j=0; j<NumImportIDs; j++) {
1040 ptr = (
double *) Imports + j*thisSizeOfPacket;
1041 jj = ToFirstPointInElementList[ImportLIDs[j]];
1042 int ElementSize = ToElementSizeList[ImportLIDs[j]];
1043 for (
int i=0; i<numVectors; i++) {
1044 if (CombineMode==
InsertAdd)
for (k=0; k<ElementSize; k++) To[i][jj+k] = 0.0;
1045 if (CombineMode==
Add || CombineMode==
InsertAdd)
for (k=0; k<ElementSize; k++) To[i][jj+k] += *ptr++;
1046 else if (CombineMode==
Insert)
for (k=0; k<ElementSize; k++) To[i][jj+k] = *ptr++;
1047 else if (CombineMode==
AbsMax) {
for (k=0; k<ElementSize; k++) { To[i][jj+k] =
EPETRA_MAX( To[i][jj+k], std::abs(*ptr)); ptr++; }}
1048 else if (CombineMode==
Epetra_Max) {
for (k=0; k<ElementSize; k++) { To[i][jj+k] =
EPETRA_MAX( To[i][jj+k], *ptr); ptr++; }}
1049 else if (CombineMode==
Epetra_Min) {
for (k=0; k<ElementSize; k++) { To[i][jj+k] =
EPETRA_MIN( To[i][jj+k], *ptr); ptr++; }}
1050 else if (CombineMode==
Average) {
for (k=0; k<ElementSize; k++) { To[i][jj+k] += *ptr++; To[i][jj+k] *= 0.5;}}
1132 double **A_Pointers =
A.Pointers();
1134 #ifdef EPETRA_HAVE_OMP_NONASSOCIATIVE 1137 const double *
const from =
Pointers_[i];
1138 const double *
const fromA = A_Pointers[i];
1140 #pragma omp parallel for reduction (+:sum) default(none) 1141 for (
int j=0; j < myLength; j++) sum += from[j] * fromA[j];
1166 double **A_Pointers =
A.Pointers();
1170 const double *
const from = A_Pointers[i];
1171 #ifdef EPETRA_HAVE_OMP 1172 #pragma omp parallel for default(none) 1174 for (
int j=0; j < myLength; j++) to[j] = std::abs(from[j]);
1185 #ifndef EPETRA_HAVE_OMP 1192 double **A_Pointers =
A.Pointers();
1196 const double *
const from = A_Pointers[i];
1197 #ifdef EPETRA_HAVE_OMP 1198 #pragma omp parallel default(none) shared(ierr) 1203 for (
int j=0; j < myLength; j++) {
1204 double value = from[j];
1207 if (value==0.0) localierr = 1;
1208 else if (localierr!=1) localierr = 2;
1214 #ifdef EPETRA_HAVE_OMP 1215 #pragma omp critical 1218 if (localierr==1) ierr = 1;
1219 else if (localierr==2 && ierr!=1) ierr = 2;
1221 #ifdef EPETRA_HAVE_OMP 1235 #ifdef EPETRA_HAVE_OMP 1238 #pragma omp parallel for default(none) shared(ScalarValue) 1239 for (
int j = 0; j < myLength; j++) to[j] = ScalarValue * to[j];
1260 double **A_Pointers = (
double**)
A.Pointers();
1264 const double *
const from = A_Pointers[i];
1265 #ifdef EPETRA_HAVE_OMP 1266 #pragma omp parallel for default(none) shared(ScalarA) 1268 for (
int j = 0; j < myLength; j++) to[j] = ScalarA * from[j];
1285 double **A_Pointers = (
double**)
A.Pointers();
1287 if (ScalarThis==0.0)
1291 const double *
const from = A_Pointers[i];
1292 #ifdef EPETRA_HAVE_OMP 1293 #pragma omp parallel for default(none) shared(ScalarA) 1295 for (
int j = 0; j < myLength; j++) to[j] = ScalarA * from[j];
1299 else if (ScalarThis==1.0)
1303 const double *
const from = A_Pointers[i];
1304 #ifdef EPETRA_HAVE_OMP 1305 #pragma omp parallel for default(none) shared(ScalarA) 1307 for (
int j = 0; j < myLength; j++) to[j] = to[j] + ScalarA * from[j];
1311 else if (ScalarA==1.0)
1315 const double *
const from = A_Pointers[i];
1316 #ifdef EPETRA_HAVE_OMP 1317 #pragma omp parallel for default(none) shared(ScalarThis) 1319 for (
int j = 0; j < myLength; j++) to[j] = ScalarThis * to[j] + from[j];
1327 const double *
const from = A_Pointers[i];
1328 #ifdef EPETRA_HAVE_OMP 1329 #pragma omp parallel for default(none) shared(ScalarA,ScalarThis) 1331 for (
int j = 0; j < myLength; j++) to[j] = ScalarThis * to[j] +
1359 if (myLength !=
A.MyLength() || myLength !=
B.MyLength())
EPETRA_CHK_ERR(-2);
1361 double **A_Pointers = (
double**)
A.Pointers();
1362 double **B_Pointers = (
double**)
B.Pointers();
1364 if (ScalarThis==0.0)
1370 const double *
const fromA = A_Pointers[i];
1371 const double *
const fromB = B_Pointers[i];
1372 #ifdef EPETRA_HAVE_OMP 1373 #pragma omp parallel for default(none) shared(ScalarB) 1375 for (
int j = 0; j < myLength; j++) to[j] = fromA[j] +
1380 else if (ScalarB==1.0)
1384 const double *
const fromA = A_Pointers[i];
1385 const double *
const fromB = B_Pointers[i];
1386 #ifdef EPETRA_HAVE_OMP 1387 #pragma omp parallel for default(none) shared(ScalarA) 1389 for (
int j = 0; j < myLength; j++) to[j] = ScalarA * fromA[j] +
1398 const double *
const fromA = A_Pointers[i];
1399 const double *
const fromB = B_Pointers[i];
1400 #ifdef EPETRA_HAVE_OMP 1401 #pragma omp parallel for default(none) shared(ScalarA,ScalarB) 1403 for (
int j = 0; j < myLength; j++) to[j] = ScalarA * fromA[j] +
1409 else if (ScalarThis==1.0)
1415 const double *
const fromA = A_Pointers[i];
1416 const double *
const fromB = B_Pointers[i];
1417 #ifdef EPETRA_HAVE_OMP 1418 #pragma omp parallel for default(none) shared(ScalarB) 1420 for (
int j = 0; j < myLength; j++) to[j] += fromA[j] +
1425 else if (ScalarB==1.0)
1429 const double *
const fromA = A_Pointers[i];
1430 const double *
const fromB = B_Pointers[i];
1431 #ifdef EPETRA_HAVE_OMP 1432 #pragma omp parallel for default(none) shared(ScalarA) 1434 for (
int j = 0; j < myLength; j++) to[j] += ScalarA * fromA[j] +
1443 const double *
const fromA = A_Pointers[i];
1444 const double *
const fromB = B_Pointers[i];
1445 #ifdef EPETRA_HAVE_OMP 1446 #pragma omp parallel for default(none) shared(ScalarA,ScalarB) 1448 for (
int j = 0; j < myLength; j++) to[j] += ScalarA * fromA[j] +
1460 const double *
const fromA = A_Pointers[i];
1461 const double *
const fromB = B_Pointers[i];
1462 #ifdef EPETRA_HAVE_OMP 1463 #pragma omp parallel for default(none) shared(ScalarA,ScalarB,ScalarThis) 1465 for (
int j = 0; j < myLength; j++) to[j] = ScalarThis * to[j] +
1471 else if (ScalarB==1.0)
1475 const double *
const fromA = A_Pointers[i];
1476 const double *
const fromB = B_Pointers[i];
1477 #ifdef EPETRA_HAVE_OMP 1478 #pragma omp parallel for default(none) shared(ScalarA,ScalarThis) 1480 for (
int j = 0; j < myLength; j++) to[j] = ScalarThis * to[j] +
1481 ScalarA * fromA[j] +
1490 const double *
const fromA = A_Pointers[i];
1491 const double *
const fromB = B_Pointers[i];
1492 #ifdef EPETRA_HAVE_OMP 1493 #pragma omp parallel for default(none) shared(ScalarA,ScalarB,ScalarThis) 1495 for (
int j = 0; j < myLength; j++) to[j] = ScalarThis * to[j] +
1496 ScalarA * fromA[j] +
1516 #ifdef EPETRA_HAVE_OMP_NONASSOCIATIVE 1519 const double *
const from =
Pointers_[i];
1521 #pragma omp parallel default(none) shared(asum) 1523 double localasum = 0.0;
1525 for (
int j=0; j< myLength; j++) localasum += std::abs(from[j]);
1526 #pragma omp critical 1559 const double *
const from =
Pointers_[i];
1561 #ifdef EPETRA_HAVE_OMP_NONASSOCIATIVE 1562 #pragma omp parallel for reduction (+:sum) default(none) 1564 for (
int j=0; j < myLength; j++) sum += from[j] * from[j];
1572 for (
int i=0; i <
NumVectors_; i++) Result[i] = std::sqrt(Result[i]);
1591 double normval = 0.0;
1592 const double *
const from =
Pointers_[i];
1593 if (myLength>0) normval = from[0];
1594 #ifdef EPETRA_HAVE_OMP 1595 #pragma omp parallel default(none) shared(normval) 1597 double localnormval = 0.0;
1599 for (
int j=0; j< myLength; j++) {
1600 localnormval =
EPETRA_MAX(localnormval,std::abs(from[j]));
1602 #pragma omp critical 1640 double *W = Weights.
Values();
1641 double **W_Pointers = Weights.
Pointers();
1645 if (!OneW) W = W_Pointers[i];
1647 const double *
const from =
Pointers_[i];
1648 #ifdef EPETRA_HAVE_OMP_NONASSOCIATIVE 1649 #pragma omp parallel for reduction (+:sum) default(none) shared(W) 1651 for (
int j=0; j < myLength; j++) {
1652 double tmp = from[j]/W[j];
1664 OneOverN = 1.0 / (double) myLength;
1666 for (
int i=0; i <
NumVectors_; i++) Result[i] = std::sqrt(Result[i]*OneOverN);
1685 const double *
const from =
Pointers_[i];
1687 if (myLength>0) MinVal = from[0];
1688 #ifdef EPETRA_HAVE_OMP 1689 #pragma omp parallel default(none) shared(MinVal) 1691 double localMinVal = MinVal;
1693 for (
int j=0; j< myLength; j++) localMinVal =
EPETRA_MIN(localMinVal,from[j]);
1694 #pragma omp critical 1701 for (
int j=0; j< myLength; j++) MinVal =
EPETRA_MIN(MinVal,from[j]);
1734 if (!epetrampicomm) {
1738 MPI_Comm mpicomm = epetrampicomm->
GetMpiComm();
1739 int numProcs = epetrampicomm->
NumProc();
1740 double* dwork =
new double[numProcs*(
NumVectors_+1)];
1749 bool overwrite = myLength == 0 ? true :
false;
1751 int myPID = epetrampicomm->
MyPID();
1752 double* dwork_ptr = dwork;
1754 for(
int j=0; j<numProcs; ++j) {
1764 double val = dwork_ptr[i];
1768 if (overwrite || (Result[i] > val)) Result[i] = val;
1773 if (overwrite) overwrite =
false;
1798 const double *
const from =
Pointers_[i];
1800 if (myLength>0) MaxVal = from[0];
1801 #ifdef EPETRA_HAVE_OMP 1802 #pragma omp parallel default(none) shared(MaxVal) 1804 double localMaxVal = MaxVal;
1806 for (
int j=0; j< myLength; j++) localMaxVal =
EPETRA_MAX(localMaxVal,from[j]);
1807 #pragma omp critical 1814 for (
int j=0; j< myLength; j++) MaxVal =
EPETRA_MAX(MaxVal,from[j]);
1847 if (!epetrampicomm) {
1851 MPI_Comm mpicomm = epetrampicomm->
GetMpiComm();
1852 int numProcs = epetrampicomm->
NumProc();
1853 double* dwork =
new double[numProcs*(
NumVectors_+1)];
1862 bool overwrite = myLength == 0 ? true :
false;
1864 int myPID = epetrampicomm->
MyPID();
1865 double* dwork_ptr = dwork;
1867 for(
int j=0; j<numProcs; ++j) {
1877 double val = dwork_ptr[i];
1881 if (overwrite || (Result[i] < val)) Result[i] = val;
1886 if (overwrite) overwrite =
false;
1916 const double *
const from =
Pointers_[i];
1917 #ifdef EPETRA_HAVE_OMP_NONASSOCIATIVE 1918 #pragma omp parallel for reduction (+:sum) default(none) 1920 for (
int j=0; j < myLength; j++) sum += from[j];
1928 for (
int i=0; i <
NumVectors_; i++) Result[i] = Result[i]*fGlobalLength;
1939 double ScalarThis ) {
1974 int A_nrows = (TransA==
'T') ?
A.NumVectors() :
A.MyLength();
1975 int A_ncols = (TransA==
'T') ?
A.MyLength() :
A.NumVectors();
1976 int B_nrows = (TransB==
'T') ?
B.NumVectors() :
B.MyLength();
1977 int B_ncols = (TransB==
'T') ?
B.MyLength() :
B.NumVectors();
1979 double Scalar_local = ScalarThis;
1982 if( myLength != A_nrows ||
1983 A_ncols != B_nrows ||
1987 bool A_is_local = (!
A.DistributedGlobal());
1988 bool B_is_local = (!
B.DistributedGlobal());
1990 bool Case1 = ( A_is_local && B_is_local && C_is_local);
1991 bool Case2 = (!A_is_local && !B_is_local && C_is_local && TransA==
'T' );
1992 bool Case3 = (!A_is_local && B_is_local && !C_is_local && TransA==
'N');
1994 if (Case2 && (!
A.Map().UniqueGIDs() || !
B.Map().UniqueGIDs())) {
EPETRA_CHK_ERR(-4);}
1999 if (Case1 || Case2 || Case3)
2001 if (ScalarThis!=0.0 && Case2)
2004 if (
MyPID!=0) Scalar_local = 0.0;
2026 double *Ap = A_tmp->
Values();
2027 double *Bp = B_tmp->
Values();
2028 double *Cp = C_tmp->
Values();
2030 GEMM(TransA, TransB, m,
n, k, ScalarAB,
2031 Ap, lda, Bp, ldb, Scalar_local, Cp, ldc);
2063 if (!
A.ConstantStride())
delete A_tmp;
2064 if (!
B.ConstantStride())
delete B_tmp;
2088 double ScalarThis) {
2094 if (ScalarAB==0.0) {
2100 if (
A.NumVectors() != 1 &&
A.NumVectors() !=
B.NumVectors())
EPETRA_CHK_ERR(-1);
2102 if (myLength !=
A.MyLength() || myLength !=
B.MyLength())
EPETRA_CHK_ERR(-3);
2104 double **A_Pointers = (
double**)
A.Pointers();
2105 double **B_Pointers = (
double**)
B.Pointers();
2108 if (
A.NumVectors() == 1 ) IncA = 0;
2110 if (ScalarThis==0.0) {
2114 const double *
const Aptr = A_Pointers[i*IncA];
2115 const double *
const Bptr = B_Pointers[i];
2117 #ifdef EPETRA_HAVE_OMP 2118 #pragma omp parallel for default(none) 2120 for (
int j = 0; j < myLength; j++) {
2121 to[j] = Aptr[j] * Bptr[j];
2129 const double *
const Aptr = A_Pointers[i*IncA];
2130 const double *
const Bptr = B_Pointers[i];
2132 #ifdef EPETRA_HAVE_OMP 2133 #pragma omp parallel for default(none) shared(ScalarAB) 2135 for (
int j = 0; j < myLength; j++) {
2136 to[j] = ScalarAB * Aptr[j] * Bptr[j];
2142 else if (ScalarThis==1.0) {
2146 const double *
const Aptr = A_Pointers[i*IncA];
2147 const double *
const Bptr = B_Pointers[i];
2149 #ifdef EPETRA_HAVE_OMP 2150 #pragma omp parallel for default(none) 2152 for (
int j = 0; j < myLength; j++) {
2153 to[j] += Aptr[j] * Bptr[j];
2160 const double *
const Aptr = A_Pointers[i*IncA];
2161 const double *
const Bptr = B_Pointers[i];
2163 #ifdef EPETRA_HAVE_OMP 2164 #pragma omp parallel for default(none) shared(ScalarAB) 2166 for (
int j = 0; j < myLength; j++) {
2167 to[j] += ScalarAB * Aptr[j] * Bptr[j];
2177 const double *
const Aptr = A_Pointers[i*IncA];
2178 const double *
const Bptr = B_Pointers[i];
2180 #ifdef EPETRA_HAVE_OMP 2181 #pragma omp parallel for default(none) shared(ScalarThis) 2183 for (
int j = 0; j < myLength; j++) {
2184 to[j] = ScalarThis * to[j] +
2193 const double *
const Aptr = A_Pointers[i*IncA];
2194 const double *
const Bptr = B_Pointers[i];
2196 #ifdef EPETRA_HAVE_OMP 2197 #pragma omp parallel for default(none) shared(ScalarThis,ScalarAB) 2199 for (
int j = 0; j < myLength; j++) {
2200 to[j] = ScalarThis * to[j] +
2201 ScalarAB * Aptr[j] * Bptr[j];
2211 double ScalarThis) {
2217 if (ScalarAB==0.0) {
2223 if (
A.NumVectors() != 1 &&
A.NumVectors() !=
B.NumVectors())
EPETRA_CHK_ERR(-1);
2225 if (myLength !=
A.MyLength() || myLength !=
B.MyLength())
EPETRA_CHK_ERR(-3);
2227 double **A_Pointers = (
double**)
A.Pointers();
2228 double **B_Pointers = (
double**)
B.Pointers();
2231 if (
A.NumVectors() == 1 ) IncA = 0;
2233 if (ScalarThis==0.0) {
2237 const double *
const Aptr = A_Pointers[i*IncA];
2238 const double *
const Bptr = B_Pointers[i];
2240 #ifdef EPETRA_HAVE_OMP 2241 #pragma omp parallel for default(none) 2243 for (
int j = 0; j < myLength; j++) {
2244 to[j] = Bptr[j] / Aptr[j];
2252 const double *
const Aptr = A_Pointers[i*IncA];
2253 const double *
const Bptr = B_Pointers[i];
2255 #ifdef EPETRA_HAVE_OMP 2256 #pragma omp parallel for default(none) shared(ScalarAB) 2258 for (
int j = 0; j < myLength; j++) {
2259 to[j] = ScalarAB * Bptr[j] / Aptr[j];
2265 else if (ScalarThis==1.0) {
2269 const double *
const Aptr = A_Pointers[i*IncA];
2270 const double *
const Bptr = B_Pointers[i];
2272 #ifdef EPETRA_HAVE_OMP 2273 #pragma omp parallel for default(none) 2275 for (
int j = 0; j < myLength; j++) {
2276 to[j] += Bptr[j] / Aptr[j];
2284 const double *
const Aptr = A_Pointers[i*IncA];
2285 const double *
const Bptr = B_Pointers[i];
2287 #ifdef EPETRA_HAVE_OMP 2288 #pragma omp parallel for default(none) shared(ScalarAB) 2290 for (
int j = 0; j < myLength; j++) {
2291 to[j] += ScalarAB * Bptr[j] / Aptr[j];
2301 const double *
const Aptr = A_Pointers[i*IncA];
2302 const double *
const Bptr = B_Pointers[i];
2304 #ifdef EPETRA_HAVE_OMP 2305 #pragma omp parallel for default(none) shared(ScalarThis) 2307 for (
int j = 0; j < myLength; j++) {
2308 to[j] = ScalarThis * to[j] +
2316 const double *
const Aptr = A_Pointers[i*IncA];
2317 const double *
const Bptr = B_Pointers[i];
2319 #ifdef EPETRA_HAVE_OMP 2320 #pragma omp parallel for default(none) shared(ScalarAB,ScalarThis) 2322 for (
int j = 0; j < myLength; j++) {
2323 to[j] = ScalarThis * to[j] + ScalarAB *
2371 if (
this != &Source)
Assign(Source);
2382 +
". The A MultiVector has NumVectors = " +
toString(
A.NumVectors()), -3);
2383 if (myLength !=
A.MyLength())
2384 throw ReportError(
"Length of MultiVectors incompatible in Assign. The this MultiVector has MyLength = " +
toString(myLength)
2385 +
". The A MultiVector has MyLength = " +
toString(
A.MyLength()), -4);
2387 double ** A_Pointers =
A.Pointers();
2390 const double *
const from = A_Pointers[i];
2391 #ifdef EPETRA_HAVE_OMP 2392 #pragma omp parallel for default(none) 2394 for (
int j=0; j<myLength; j++) to[j] = from[j];
2404 double * source = 0;
2405 if (myLength>0) source =
new double[myLength*
NumVectors_];
2406 double * target = 0;
2413 double * tmp1 = source;
2416 for (
int j=0; j< myLength; j++) *tmp1++ = *tmp2++;
2418 if (myLength>0) target =
new double[myLength*
NumVectors_];
2422 if (myLength>0)
delete [] source;
2424 double * tmp2 = target;
2427 for (
int j=0; j< myLength; j++) *tmp1++ = *tmp2++;
2429 if (myLength>0)
delete [] target;
2451 for (
int iproc=0; iproc <
NumProc; iproc++) {
2454 int NumMyElements1 =
Map(). NumMyElements();
2456 int * FirstPointInElementList1 = NULL;
2462 os <<
" MyPID"; os <<
" ";
2464 if (MaxElementSize1==1)
2468 for (
int j = 0; j < NumVectors1 ; j++)
2475 for (
int i=0; i < NumMyElements1; i++) {
2479 os <<
MyPID; os <<
" ";
2481 if (MaxElementSize1==1) {
2482 if(
Map().GlobalIndicesInt())
2484 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 2486 os << MyGlobalElements1[i] <<
" ";
2488 throw ReportError(
"Epetra_MultiVector::Print: ERROR, GlobalIndicesInt but no API for it.",-1);
2491 else if(
Map().GlobalIndicesLongLong())
2493 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 2495 os << MyGlobalElements1[i] <<
" ";
2497 throw ReportError(
"Epetra_MultiVector::Print: ERROR, GlobalIndicesLongLong but no API for it.",-1);
2501 throw ReportError(
"Epetra_MultiVector::Print ERROR, Don't know map global index type.",-1);
2506 if(
Map().GlobalIndicesInt())
2508 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 2510 os << MyGlobalElements1[i]<<
"/" << ii <<
" ";
2512 throw ReportError(
"Epetra_MultiVector::Print: ERROR, GlobalIndicesInt but no API for it.",-1);
2515 else if(
Map().GlobalIndicesLongLong())
2517 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 2519 os << MyGlobalElements1[i]<<
"/" << ii <<
" ";
2521 throw ReportError(
"Epetra_MultiVector::Print: ERROR, GlobalIndicesLongLong but no API for it.",-1);
2525 throw ReportError(
"Epetra_MultiVector::Print ERROR, Don't know map global index type.",-1);
2527 iii = FirstPointInElementList1[i]+ii;
2529 for (
int j = 0; j < NumVectors1 ; j++)
2532 os << A_Pointers[j][iii];
int ChangeMyValue(int MyBlockRow, int BlockRowOffset, int VectorIndex, double ScalarValue, bool SumInto)
int Norm2(double *Result) const
Compute 2-norm of each vector in multi-vector.
int Norm1(double *Result) const
Compute 1-norm of each vector in multi-vector.
Epetra_MultiVector: A class for constructing and using dense multi-vectors, vectors and matrices in p...
double * Values() const
Get pointer to MultiVector values.
int NumProc() const
Returns total number of processes.
const Epetra_BlockMap & Map() const
Returns the address of the Epetra_BlockMap for this multi-vector.
int MyGlobalElements(int *MyGlobalElementList) const
Puts list of global elements on this processor into the user-provided array.
int Abs(const Epetra_MultiVector &A)
Puts element-wise absolute values of input Multi-vector in target.
int Dot(const Epetra_MultiVector &A, double *Result) const
Computes dot product of each corresponding pair of vectors.
int Random()
Set multi-vector values to random numbers.
int ReplaceMyValue(int MyRow, int VectorIndex, double ScalarValue)
Replace current value at the specified (MyRow, VectorIndex) location with ScalarValue.
void UpdateFlops(int Flops_in) const
Increment Flop count for this object.
bool UniqueGIDs() const
Returns true if map GIDs are 1-to-1.
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.
int MyLength() const
Returns the local vector length on the calling processor of vectors in the multi-vector.
Epetra_Distributor: The Epetra Gather/Scatter Setup Base Class.
int AllocateForCopy(void)
int MeanValue(double *Result) const
Compute mean (average) value of each vector in multi-vector.
virtual void Print(std::ostream &os) const
Print method.
bool ConstantElementSize() const
Returns true if map has constant element size.
int NormWeighted(const Epetra_MultiVector &Weights, double *Result) const
Compute Weighted 2-norm (RMS Norm) of each vector in multi-vector.
float DOT(const int N, const float *X, const float *Y, const int INCX=1, const int INCY=1) const
Epetra_BLAS dot product function (SDOT).
Epetra_OffsetIndex: This class builds index for efficient mapping of data from one Epetra_CrsGraph ba...
int IAMAX(const int N, const float *X, const int INCX=1) const
Epetra_BLAS arg maximum of absolute value function (ISAMAX)
float ASUM(const int N, const float *X, const int INCX=1) const
Epetra_BLAS one norm function (SASUM).
int PutScalar(double ScalarConstant)
Initialize all values in a multi-vector with constant value.
int MyPID() const
Return my process ID.
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().
#define EPETRA_CHK_ERR(a)
int NormInf(double *Result) const
Compute Inf-norm of each vector in multi-vector.
int Stride() const
Returns the stride between vectors in the multi-vector (only meaningful if ConstantStride() is true)...
void UpdateDoubleTemp() const
int MaxValue(double *Result) const
Compute maximum value of each vector in multi-vector.
int * FirstPointInElementList() const
Pointer to internal array containing a mapping between the local elements and the first local point n...
Epetra_Vector: A class for constructing and using dense vectors on a parallel computer.
int ExtractCopy(double *A, int MyLDA) const
Put multi-vector values into user-provided two-dimensional array.
virtual int ReportError(const std::string Message, int ErrorCode) const
Error reporting method.
Epetra_MultiVector(const Epetra_BlockMap &Map, int NumVectors, bool zeroOut=true)
Basic Epetra_MultiVector constuctor.
Epetra_MpiComm: The Epetra MPI Communication Class.
int FirstPointInElement(int LID) const
Returns the requested entry in the FirstPointInElementList; see FirstPointInElementList() for details...
virtual void Barrier() const =0
Epetra_Comm Barrier function.
int MinValue(double *Result) const
Compute minimum value of each vector in multi-vector.
int NumProc() const
Returns total number of processes.
long long * MyGlobalElements64() const
virtual int MyPID() const =0
Return my process ID.
int * ElementSizeList() const
List of the element sizes corresponding to the array MyGlobalElements().
virtual int MaxAll(double *PartialMaxs, double *GlobalMaxs, int Count) const =0
Epetra_Comm Global Max function.
Epetra_Vector ** Vectors_
int SumIntoGlobalValue(int GlobalRow, int VectorIndex, double ScalarValue)
Adds ScalarValue to existing value at the specified (GlobalRow, VectorIndex) location.
virtual int SumAll(double *PartialSums, double *GlobalSums, int Count) const =0
Epetra_Comm Global Sum function.
Epetra_MultiVector & operator=(const Epetra_MultiVector &Source)
= Operator.
int Scale(double ScalarValue)
Scale the current values of a multi-vector, this = ScalarValue*this.
virtual ~Epetra_MultiVector()
Epetra_MultiVector destructor.
int SetSeed(unsigned int Seed_in)
Set seed for Random function.
const Epetra_Comm * Comm_
int ResetView(double **ArrayOfPointers)
Reset the view of an existing multivector to point to new user data.
int NumVectors() const
Returns the number of vectors in the multi-vector.
int NumMyElements() const
Number of elements on the calling processor.
const double Epetra_MaxDouble
int ReplaceGlobalValue(int GlobalRow, int VectorIndex, double ScalarValue)
Replace current value at the specified (GlobalRow, VectorIndex) location with ScalarValue.
Epetra_CompObject: Functionality and data that is common to all computational classes.
MPI_Comm GetMpiComm() const
Get the MPI Communicator (identical to Comm() method; used when we know we are MPI.
int CheckSizes(const Epetra_SrcDistObject &A)
Allows the source and target (this) objects to be compared for compatibility, return nonzero if not...
int ElementSize() const
Returns the size of elements in the map; only valid if map has constant element size.
std::string toString(const int &x) const
int SumIntoMyValue(int MyRow, int VectorIndex, double ScalarValue)
Adds ScalarValue to existing value at the specified (MyRow, VectorIndex) location.
Epetra_BlockMap: A class for partitioning block element vectors and matrices.
const Epetra_Comm & Comm() const
Access function for Epetra_Comm communicator.
void GEMM(const char TRANSA, const char TRANSB, const int M, const int N, const int K, const float ALPHA, const float *A, const int LDA, const float *B, const int LDB, const float BETA, float *C, const int LDC) const
Epetra_BLAS matrix-matrix multiply function (SGEMM)
int MyPID() const
Return my process ID.
void UpdateVectors() const
const double Epetra_MinDouble
int Update(double ScalarA, const Epetra_MultiVector &A, double ScalarThis)
Update multi-vector values with scaled values of A, this = ScalarThis*this + ScalarA*A.
int Multiply(char TransA, char TransB, double ScalarAB, const Epetra_MultiVector &A, const Epetra_MultiVector &B, double ScalarThis)
Matrix-Matrix multiplication, this = ScalarThis*this + ScalarAB*A*B.
virtual int NumProc() const =0
Returns total number of processes.
double RandomDouble()
Returns a random double on the interval (-1.0,1.0)
int ChangeGlobalValue(int_type GlobalBlockRow, int BlockRowOffset, int VectorIndex, double ScalarValue, bool SumInto)
void Assign(const Epetra_MultiVector &rhs)
int AllocateForView(void)
int MaxElementSize() const
Maximum element size across all processors.
int Reciprocal(const Epetra_MultiVector &A)
Puts element-wise reciprocal values of input Multi-vector in target.
int ExtractView(double **A, int *MyLDA) const
Set user-provided addresses of A and MyLDA.
int ReplaceMap(const Epetra_BlockMap &map)
Replace map, only if new map has same point-structure as current map.
Epetra_SrcDistObject: A class for supporting flexible source distributed objects for import/export op...
int NumMyPoints() const
Number of local points for this map; equals the sum of all element sizes on the calling processor...
void SCAL(const int N, const float ALPHA, float *X, const int INCX=1) const
Epetra_BLAS vector scale function (SSCAL)
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().
double ** Pointers() const
Get pointer to individual vector pointers.
bool DistributedGlobal() const
Returns true if this multi-vector is distributed global, i.e., not local replicated.
Epetra_DistObject: A class for constructing and using dense multi-vectors, vectors and matrices in pa...
Epetra_Vector *& operator()(int i)
Vector access function.
int ReciprocalMultiply(double ScalarAB, const Epetra_MultiVector &A, const Epetra_MultiVector &B, double ScalarThis)
Multiply a Epetra_MultiVector by the reciprocal of another, element-by-element.