43 #include "Epetra_ConfigDefs.h"
44 #include "Epetra_Map.h"
45 #include "Epetra_Util.h"
46 #include "Epetra_Export.h"
47 #include "Epetra_Import.h"
48 #include "Epetra_MultiVector.h"
49 #include "Epetra_Vector.h"
50 #include "Epetra_GIDTypeVector.h"
51 #include "Epetra_Comm.h"
52 #include "Epetra_LinearProblem.h"
53 #include "Epetra_MapColoring.h"
89 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
125 std::cout <<
"\nAnalyzed Singleton Problem:\n";
126 std::cout <<
"---------------------------\n";
130 std::cout <<
"Singletons Detected!" << std::endl;;
131 std::cout <<
"Num Singletons: " <<
NumSingletons() << std::endl;
136 std::cout <<
"No Singletons Detected!" << std::endl;
151 std::cout <<
"---------------------------\n\n";
176 std::cout <<
"\nConstructedSingleton Problem:\n";
177 std::cout <<
"---------------------------\n";
180 std::cout <<
"---------------------------\n\n";
190 if( ierr ) std::cout <<
"EDT_LinearProblem_CrsSingletonFilter::UpdateReducedProblem FAILED!\n";
199 if( ierr ) std::cout <<
"EDT_LinearProblem_CrsSingletonFilter::ComputeFullSolution FAILED!\n";
247 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
263 if (fullMatrix==0) EPETRA_CHK_ERR(-2);
264 if (fullMatrix->NumGlobalRows64()==0) EPETRA_CHK_ERR(-3);
265 if (fullMatrix->NumGlobalNonzeros64()==0) EPETRA_CHK_ERR(-4);
291 for (
int i=0; i<NumMyRows; i++) {
293 EPETRA_CHK_ERR(
GetRow(i, NumIndices, Indices));
294 for (
int j=0; j<NumIndices; j++) {
295 int ColumnIndex = Indices[j];
296 ColProfiles[ColumnIndex]++;
299 RowIDs[ColumnIndex] = i;
304 ColHasRowWithSingleton[j2]++;
306 colMapColors[j2] = 1;
339 if (ColHasRowWithSingleton.
MaxValue()>1) {
348 for (
int j=0; j<NumMyCols; j++) {
351 if (ColProfiles[j]==1 ) {
355 if (rowMapColors[i2]!=1) {
356 RowHasColWithSingleton[i2]++;
357 rowMapColors[i2] = 2;
363 EPETRA_CHK_ERR(
GetRow(i2, NumIndices, Indices));
364 for (
int jj=0; jj<NumIndices; jj++) {
365 NewColProfiles[Indices[jj]]--;
370 else if (ColHasRowWithSingleton[j]==1 && rowMapColors[i2]!=1) {
374 if (RowHasColWithSingleton.MaxValue()>1) {
380 ColHasRowWithSingleton));
382 for (
int i=0; i<NumMyRows; i++) {
383 if (rowMapColors[i]==2) rowMapColors[i] = 1;
392 template<
typename int_type>
393 int LinearProblem_CrsSingletonFilter::TConstructReducedProblem(
Epetra_LinearProblem * Problem) {
396 if (Problem==0) EPETRA_CHK_ERR(-2);
401 if (Problem->
GetRHS()==0) EPETRA_CHK_ERR(-4);
402 if (Problem->
GetLHS()==0) EPETRA_CHK_ERR(-5);
445 OrigReducedMatrixDomainMap_ = 0;
452 int NumVectors = FullLHS->NumVectors();
467 int ColSingletonCounter = 0;
468 for (i=0; i<NumMyRows; i++) {
472 EPETRA_CHK_ERR(
GetRowGCIDs(i, NumEntries, Values, Indices));
480 if (ierr<0) EPETRA_CHK_ERR(ierr);
484 EPETRA_CHK_ERR(
GetRow(i, NumEntries, Values, Indices));
486 double pivot = Values[0];
487 if (pivot==0.0) EPETRA_CHK_ERR(-1);
488 int indX = Indices[0];
489 for (j=0; j<NumVectors; j++)
496 for (j=0; j<NumEntries; j++) {
497 if (Indices[j]==targetCol) {
498 double pivot = Values[j];
499 if (pivot==0.0) EPETRA_CHK_ERR(-2);
502 ColSingletonCounter++;
520 FullLHS->PutScalar(0.0);
542 EPETRA_CHK_ERR(
tempB_->Update(1.0, *FullRHS, -1.0));
553 ReducedProblem_ = Teuchos::rcp( Problem,
false );
554 ReducedMatrix_ = Teuchos::rcp( dynamic_cast<Epetra_CrsMatrix *>(Problem->
GetMatrix()),
false );
557 double fn = (double)
FullMatrix()->NumGlobalRows64();
558 double fnnz = (double)
FullMatrix()->NumGlobalNonzeros64();
560 double rnnz = (double)
ReducedMatrix()->NumGlobalNonzeros64();
570 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
572 return TConstructReducedProblem<int>(Problem);
576 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
578 return TConstructReducedProblem<long long>(Problem);
582 throw "LinearProblem_CrsSingletonFilter::ConstructReducedProblem: ERROR, GlobalIndices type unknown.";
586 template<
typename int_type>
591 if (Problem==0) EPETRA_CHK_ERR(-1);
596 if (Problem->
GetRHS()==0) EPETRA_CHK_ERR(-3);
597 if (Problem->
GetLHS()==0) EPETRA_CHK_ERR(-4);
606 int NumVectors = FullLHS->NumVectors();
612 int ColSingletonCounter = 0;
613 for (i=0; i<NumMyRows; i++) {
617 EPETRA_CHK_ERR(
GetRowGCIDs(i, NumEntries, Values, Indices));
624 if (ierr<0) EPETRA_CHK_ERR(ierr);
629 EPETRA_CHK_ERR(
GetRow(i, NumEntries, Values, Indices));
631 double pivot = Values[0];
632 if (pivot==0.0) EPETRA_CHK_ERR(-1);
633 int indX = Indices[0];
634 for (j=0; j<NumVectors; j++)
641 double pivot = Values[j];
642 if (pivot==0.0) EPETRA_CHK_ERR(-2);
644 ColSingletonCounter++;
655 FullLHS->PutScalar(0.0);
677 EPETRA_CHK_ERR(
tempB_->Update(1.0, *FullRHS, -1.0));
679 ReducedRHS_->PutScalar(0.0);
685 ReducedProblem_ = Teuchos::rcp( Problem,
false );
686 ReducedMatrix_ = Teuchos::rcp( dynamic_cast<Epetra_CrsMatrix *>(Problem->
GetMatrix()),
false );
693 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
695 return TUpdateReducedProblem<int>(Problem);
699 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
701 return TUpdateReducedProblem<long long>(Problem);
705 throw "LinearProblem_CrsSingletonFilter::UpdateReducedProblem: ERROR, GlobalIndices type unknown.";
708 template<
typename int_type>
709 int LinearProblem_CrsSingletonFilter::TConstructRedistributeExporter(
Epetra_Map * SourceMap,
Epetra_Map * TargetMap,
713 int_type IndexBase = (int_type) SourceMap->IndexBase64();
714 if (IndexBase!=(int_type) TargetMap->IndexBase64()) EPETRA_CHK_ERR(-1);
722 Epetra_Map ContiguousTargetMap((int_type) -1, TargetNumMyElements, IndexBase,Comm);
725 Epetra_Map ContiguousSourceMap((int_type) -1, SourceNumMyElements, IndexBase, Comm);
727 assert(ContiguousSourceMap.NumGlobalElements64()==ContiguousTargetMap.NumGlobalElements64());
730 int_type* SourceMapMyGlobalElements = 0;
731 SourceMap->MyGlobalElementsPtr(SourceMapMyGlobalElements);
735 Epetra_Export Exporter(ContiguousSourceMap, ContiguousTargetMap);
739 TargetIndices.Export(SourceIndices, Exporter, Insert);
743 RedistributeMap =
new Epetra_Map((int_type) -1, TargetNumMyElements, TargetIndices.Values(), IndexBase, Comm);
746 RedistributeExporter =
new Epetra_Export(*SourceMap, *RedistributeMap);
753 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
755 return TConstructRedistributeExporter<int>(SourceMap, TargetMap, RedistributeExporter, RedistributeMap);
759 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
761 return TConstructRedistributeExporter<long long>(SourceMap, TargetMap, RedistributeExporter, RedistributeMap);
765 throw "LinearProblem_CrsSingletonFilter::ConstructRedistributeExporter: ERROR, GlobalIndices type unknown.";
780 FullLHS->Update(1.0, *
tempX_, 1.0);
791 int NumVectors =
tempB_->NumVectors();
796 for (jj=0; jj<NumVectors; jj++)
797 (*
tempExportX_)[jj][j]= ((*FullRHS)[jj][i] - (*tempB_)[jj][i])/pivot;
809 FullLHS->Update(1.0, *
tempX_, 1.0);
823 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
835 EPETRA_CHK_ERR(
FullCrsMatrix()->Graph().ExtractMyRowView(Row, NumIndices, Indices));
846 double * & Values,
int * & Indices) {
860 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
862 double * & Values,
int * & GlobalIndices) {
873 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
875 double * & Values,
long long * & GlobalIndices) {
907 int NumMyColSingletonstmp = 0;
908 for (j=0; j<NumMyCols; j++) {
910 if ( ColProfiles[j]==1 && rowMapColors[i]!=1 ) {
913 NumMyColSingletonstmp++;
917 else if (NewColProfiles[j]==0 && ColHasRowWithSingleton[j]!=1 && rowMapColors[i]==0) {
virtual int ReplaceGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
virtual const Epetra_Map & RowMatrixRowMap() const =0
Epetra_MultiVector * GetLHS() const
Epetra_MultiVector * GetRHS() const
bool GlobalIndicesLongLong() const
bool SameAs(const Epetra_BlockMap &Map) const
bool GlobalIndicesInt() const
Epetra_Map * GenerateMap(int Color) const
virtual int MyPID() const =0
virtual int NumMyCols() const =0
int ExtractMyRowView(int MyRow, int &NumEntries, double *&Values, int *&Indices) const
virtual int SumAll(double *PartialSums, double *GlobalSums, int Count) const =0
int NumMyElements() const
virtual int MaxNumEntries() const =0
virtual const Epetra_Comm & Comm() const =0
virtual const Epetra_Import * RowMatrixImporter() const =0
virtual int NumMyRows() const =0
const Epetra_Comm & Comm() const
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
Epetra_RowMatrix * GetMatrix() const
static void EPETRA_LIB_DLL_EXPORT Sort(bool SortAscending, int NumKeys, T *Keys, int NumDoubleCompanions, double **DoubleCompanions, int NumIntCompanions, int **IntCompanions, int NumLongLongCompanions, long long **LongLongCompanions)