43 #include "Ifpack_ConfigDefs.h" 44 #include "Ifpack_Preconditioner.h" 46 #include "Epetra_Comm.h" 47 #include "Epetra_CrsMatrix.h" 48 #include "Epetra_CrsGraph.h" 49 #include "Epetra_Map.h" 50 #include "Epetra_BlockMap.h" 51 #include "Epetra_Import.h" 52 #include "Epetra_MultiVector.h" 53 #include "Epetra_Vector.h" 57 std::cout <<
"================================================================================" << std::endl;
68 if (Comm.MyPID() == 0) cout <<
"Host and Process Ids for tasks" << endl;
69 for (
int i = 0; i <Comm.NumProc() ; i++) {
70 if (i == Comm.MyPID() ) {
71 #if defined(TFLOP) || defined(JANUS_STLPORT) 72 sprintf(buf,
"Host: %s PID: %d",
"janus", getpid());
74 sprintf(buf,
"Windows compiler, unknown hostname and PID!");
76 gethostname(hostname,
sizeof(hostname));
77 sprintf(buf,
"Host: %s\tComm.MyPID(): %d\tPID: %d",
78 hostname, Comm.MyPID(), getpid());
82 #if !( defined(_WIN32) ) 87 if(Comm.MyPID() == 0) {
89 printf(
"** Pausing to attach debugger...\n");
90 printf(
"** You may now attach debugger to the processes listed above.\n");
92 printf(
"** Enter a character to continue > "); fflush(stdout);
94 TEUCHOS_ASSERT(scanf(
"%c",&go) != EOF);
103 const int OverlappingLevel)
106 if (OverlappingLevel == 0)
108 if (Matrix->Comm().NumProc() == 1)
111 Epetra_CrsMatrix* OverlappingMatrix;
112 OverlappingMatrix = 0;
113 Epetra_Map* OverlappingMap;
114 OverlappingMap = (Epetra_Map*)&(Matrix->RowMatrixRowMap());
116 const Epetra_RowMatrix* OldMatrix;
117 const Epetra_Map* DomainMap = &(Matrix->OperatorDomainMap());
118 const Epetra_Map* RangeMap = &(Matrix->OperatorRangeMap());
120 for (
int level = 1; level <= OverlappingLevel ; ++level) {
122 if (OverlappingMatrix)
123 OldMatrix = OverlappingMatrix;
127 Epetra_Import* OverlappingImporter;
128 OverlappingImporter = (Epetra_Import*)OldMatrix->RowMatrixImporter();
129 int NumMyElements = OverlappingImporter->TargetMap().NumMyElements();
134 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 135 if(OverlappingImporter->TargetMap().GlobalIndicesInt()) {
136 int* MyGlobalElements = OverlappingImporter->TargetMap().MyGlobalElements();
137 OverlappingMap =
new Epetra_Map(-1,NumMyElements,MyGlobalElements,
142 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 143 if(OverlappingImporter->TargetMap().GlobalIndicesLongLong()) {
144 long long* MyGlobalElements = OverlappingImporter->TargetMap().MyGlobalElements64();
145 OverlappingMap =
new Epetra_Map((
long long) -1,NumMyElements,MyGlobalElements,
150 throw "Ifpack_CreateOverlappingCrsMatrix: GlobalIndices type unknown";
152 if (level < OverlappingLevel)
153 OverlappingMatrix =
new Epetra_CrsMatrix(Copy, *OverlappingMap, 0);
158 OverlappingMatrix =
new Epetra_CrsMatrix(Copy, *OverlappingMap,
161 OverlappingMatrix->Import(*OldMatrix, *OverlappingImporter, Insert);
162 if (level < OverlappingLevel) {
163 OverlappingMatrix->FillComplete(*DomainMap, *RangeMap);
166 OverlappingMatrix->FillComplete(*DomainMap, *RangeMap);
169 delete OverlappingMap;
174 OverlappingMatrix->FillComplete();
178 return(OverlappingMatrix);
183 const int OverlappingLevel)
186 if (OverlappingLevel == 0)
188 if (Graph->Comm().NumProc() == 1)
191 Epetra_CrsGraph* OverlappingGraph;
192 Epetra_BlockMap* OverlappingMap;
193 OverlappingGraph =
const_cast<Epetra_CrsGraph*
>(Graph);
194 OverlappingMap =
const_cast<Epetra_BlockMap*
>(&(Graph->RowMap()));
196 Epetra_CrsGraph* OldGraph;
197 Epetra_BlockMap* OldMap;
198 const Epetra_BlockMap* DomainMap = &(Graph->DomainMap());
199 const Epetra_BlockMap* RangeMap = &(Graph->RangeMap());
201 for (
int level = 1; level <= OverlappingLevel ; ++level) {
203 OldGraph = OverlappingGraph;
204 OldMap = OverlappingMap;
206 Epetra_Import* OverlappingImporter;
207 OverlappingImporter =
const_cast<Epetra_Import*
>(OldGraph->Importer());
208 OverlappingMap =
new Epetra_BlockMap(OverlappingImporter->TargetMap());
210 if (level < OverlappingLevel)
211 OverlappingGraph =
new Epetra_CrsGraph(Copy, *OverlappingMap, 0);
216 OverlappingGraph =
new Epetra_CrsGraph(Copy, *OverlappingMap,
219 OverlappingGraph->Import(*OldGraph, *OverlappingImporter, Insert);
220 if (level < OverlappingLevel)
221 OverlappingGraph->FillComplete(*DomainMap, *RangeMap);
224 OverlappingImporter =
new Epetra_Import(*OverlappingMap, *DomainMap);
225 OverlappingGraph->FillComplete(*DomainMap, *RangeMap);
233 delete OverlappingMap;
234 OverlappingGraph->FillComplete();
238 return(OverlappingGraph);
246 return std::string(s);
254 return std::string(s);
259 const Epetra_MultiVector& X,
const Epetra_MultiVector&Y)
264 if (X.Comm().MyPID() == 0) {
265 cout <<
"***** " << Label << endl;
274 const Epetra_MultiVector& X,
const Epetra_MultiVector&Y)
279 Epetra_MultiVector RHS(X);
280 std::vector<double> Norm2;
281 Norm2.resize(X.NumVectors());
283 IFPACK_CHK_ERR(A.Multiply(
false,X,RHS));
284 RHS.Update(1.0, Y, -1.0);
286 RHS.Norm2(&Norm2[0]);
288 if (X.Comm().MyPID() == 0) {
289 cout <<
"***** iter: " << iter <<
": ||Ax - b||_2 = " 299 void Ifpack_PrintSparsity_Simple(
const Epetra_RowMatrix& A)
304 int MaxEntries = A.MaxNumEntries();
305 std::vector<int> Indices(MaxEntries);
306 std::vector<double> Values(MaxEntries);
307 std::vector<bool> FullRow(A.NumMyRows());
310 for (
int j = 0 ; j < A.NumMyRows() ; ++j)
312 cout <<
"-+" << endl;
314 for (
int i = 0 ; i < A.NumMyRows() ; ++i) {
317 A.ExtractMyRowCopy(i,MaxEntries,Length,
318 &Values[0], &Indices[0]);
320 for (
int j = 0 ; j < A.NumMyRows() ; ++j)
323 for (
int j = 0 ; j < Length ; ++j) {
324 FullRow[Indices[j]] =
true;
328 for (
int j = 0 ; j < A.NumMyRows() ; ++j) {
334 cout <<
" |" << endl;
338 for (
int j = 0 ; j < A.NumMyRows() ; ++j)
340 cout <<
"-+" << endl << endl;
346 double Ifpack_FrobeniusNorm(
const Epetra_RowMatrix& A)
348 double MyNorm = 0.0, GlobalNorm;
350 std::vector<int> colInd(A.MaxNumEntries());
351 std::vector<double> colVal(A.MaxNumEntries());
353 for (
int i = 0 ; i < A.NumMyRows() ; ++i) {
356 IFPACK_CHK_ERR(A.ExtractMyRowCopy(i,A.MaxNumEntries(),Nnz,
357 &colVal[0],&colInd[0]));
359 for (
int j = 0 ; j < Nnz ; ++j) {
360 MyNorm += colVal[j] * colVal[j];
364 A.Comm().SumAll(&MyNorm,&GlobalNorm,1);
366 return(sqrt(GlobalNorm));
376 static void print(
const char str[], T val)
378 std::cout.width(30); std::cout.setf(std::ios::left);
380 std::cout <<
" = " << val << std::endl;
384 static void print(
const char str[], T val,
double percentage)
386 std::cout.width(30); std::cout.setf(std::ios::left);
389 std::cout.width(20); std::cout.setf(std::ios::left);
391 std::cout <<
" ( " << percentage <<
" %)" << std::endl;
394 static void print(
const char str[], T one, T two, T three,
bool equal =
true)
398 std::cout.width(30); std::cout.setf(std::ios::left);
404 std::cout.width(15); std::cout.setf(std::ios::left);
406 std::cout.width(15); std::cout.setf(std::ios::left);
408 std::cout.width(15); std::cout.setf(std::ios::left);
416 #include "Epetra_FECrsMatrix.h" 419 const int NumPDEEqns)
422 int NumMyRows = A.NumMyRows();
423 long long NumGlobalRows = A.NumGlobalRows64();
424 long long NumGlobalCols = A.NumGlobalCols64();
425 long long MyBandwidth = 0, GlobalBandwidth;
426 long long MyLowerNonzeros = 0, MyUpperNonzeros = 0;
427 long long GlobalLowerNonzeros, GlobalUpperNonzeros;
428 long long MyDiagonallyDominant = 0, GlobalDiagonallyDominant;
429 long long MyWeaklyDiagonallyDominant = 0, GlobalWeaklyDiagonallyDominant;
430 double MyMin, MyAvg, MyMax;
431 double GlobalMin, GlobalAvg, GlobalMax;
432 long long GlobalStorage;
434 bool verbose = (A.Comm().MyPID() == 0);
436 GlobalStorage =
sizeof(
int*) * NumGlobalRows +
437 sizeof(
int) * A.NumGlobalNonzeros64() +
438 sizeof(double) * A.NumGlobalNonzeros64();
443 print<const char*>(
"Label", A.Label());
444 print<long long>(
"Global rows", NumGlobalRows);
445 print<long long>(
"Global columns", NumGlobalCols);
446 print<long long>(
"Stored nonzeros", A.NumGlobalNonzeros64());
447 print<long long>(
"Nonzeros / row", A.NumGlobalNonzeros64() / NumGlobalRows);
448 print<double>(
"Estimated storage (Mbytes)", 1.0e-6 * GlobalStorage);
451 long long NumMyActualNonzeros = 0, NumGlobalActualNonzeros;
452 long long NumMyEmptyRows = 0, NumGlobalEmptyRows;
453 long long NumMyDirichletRows = 0, NumGlobalDirichletRows;
455 std::vector<int> colInd(A.MaxNumEntries());
456 std::vector<double> colVal(A.MaxNumEntries());
458 Epetra_Vector Diag(A.RowMatrixRowMap());
459 Epetra_Vector RowSum(A.RowMatrixRowMap());
461 RowSum.PutScalar(0.0);
463 for (
int i = 0 ; i < NumMyRows ; ++i) {
465 long long GRID = A.RowMatrixRowMap().GID64(i);
467 IFPACK_CHK_ERR(A.ExtractMyRowCopy(i,A.MaxNumEntries(),Nnz,
468 &colVal[0],&colInd[0]));
474 NumMyDirichletRows++;
476 for (
int j = 0 ; j < Nnz ; ++j) {
478 double v = colVal[j];
480 if (colVal[j] != 0.0)
481 NumMyActualNonzeros++;
483 long long GCID = A.RowMatrixColMap().GID64(colInd[j]);
492 else if (GCID > GRID)
494 long long b = GCID - GRID;
500 if (Diag[i] > RowSum[i])
501 MyDiagonallyDominant++;
503 if (Diag[i] >= RowSum[i])
504 MyWeaklyDiagonallyDominant++;
506 RowSum[i] += Diag[i];
513 A.Comm().SumAll(&MyDiagonallyDominant,&GlobalDiagonallyDominant,1);
514 A.Comm().SumAll(&MyWeaklyDiagonallyDominant,&GlobalWeaklyDiagonallyDominant,1);
515 A.Comm().SumAll(&NumMyActualNonzeros, &NumGlobalActualNonzeros, 1);
516 A.Comm().SumAll(&NumMyEmptyRows, &NumGlobalEmptyRows, 1);
517 A.Comm().SumAll(&NumMyDirichletRows, &NumGlobalDirichletRows, 1);
518 A.Comm().SumAll(&MyBandwidth, &GlobalBandwidth, 1);
519 A.Comm().SumAll(&MyLowerNonzeros, &GlobalLowerNonzeros, 1);
520 A.Comm().SumAll(&MyUpperNonzeros, &GlobalUpperNonzeros, 1);
521 A.Comm().SumAll(&MyDiagonallyDominant, &GlobalDiagonallyDominant, 1);
522 A.Comm().SumAll(&MyWeaklyDiagonallyDominant, &GlobalWeaklyDiagonallyDominant, 1);
524 double NormOne = A.NormOne();
525 double NormInf = A.NormInf();
526 double NormF = Ifpack_FrobeniusNorm(A);
530 print<long long>(
"Actual nonzeros", NumGlobalActualNonzeros);
531 print<long long>(
"Nonzeros in strict lower part", GlobalLowerNonzeros);
532 print<long long>(
"Nonzeros in strict upper part", GlobalUpperNonzeros);
534 print<long long>(
"Empty rows", NumGlobalEmptyRows,
535 100.0 * NumGlobalEmptyRows / NumGlobalRows);
536 print<long long>(
"Dirichlet rows", NumGlobalDirichletRows,
537 100.0 * NumGlobalDirichletRows / NumGlobalRows);
538 print<long long>(
"Diagonally dominant rows", GlobalDiagonallyDominant,
539 100.0 * GlobalDiagonallyDominant / NumGlobalRows);
540 print<long long>(
"Weakly diag. dominant rows",
541 GlobalWeaklyDiagonallyDominant,
542 100.0 * GlobalWeaklyDiagonallyDominant / NumGlobalRows);
544 print<long long>(
"Maximum bandwidth", GlobalBandwidth);
547 print(
"",
"one-norm",
"inf-norm",
"Frobenius",
false);
548 print(
"",
"========",
"========",
"=========",
false);
551 print<double>(
"A", NormOne, NormInf, NormF);
554 if (Cheap ==
false) {
558 Epetra_FECrsMatrix AplusAT(Copy, A.RowMatrixRowMap(), 0);
559 Epetra_FECrsMatrix AminusAT(Copy, A.RowMatrixRowMap(), 0);
561 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 562 if(A.RowMatrixRowMap().GlobalIndicesInt()) {
563 for (
int i = 0 ; i < NumMyRows ; ++i) {
565 int GRID = A.RowMatrixRowMap().GID(i);
569 IFPACK_CHK_ERR(A.ExtractMyRowCopy(i,A.MaxNumEntries(),Nnz,
570 &colVal[0],&colInd[0]));
572 for (
int j = 0 ; j < Nnz ; ++j) {
574 int GCID = A.RowMatrixColMap().GID(colInd[j]);
577 double plus_val = colVal[j];
578 double minus_val = -colVal[j];
580 if (AplusAT.SumIntoGlobalValues(1,&GRID,1,&GCID,&plus_val) != 0) {
581 IFPACK_CHK_ERR(AplusAT.InsertGlobalValues(1,&GRID,1,&GCID,&plus_val));
584 if (AplusAT.SumIntoGlobalValues(1,&GCID,1,&GRID,&plus_val) != 0) {
585 IFPACK_CHK_ERR(AplusAT.InsertGlobalValues(1,&GCID,1,&GRID,&plus_val));
588 if (AminusAT.SumIntoGlobalValues(1,&GRID,1,&GCID,&plus_val) != 0) {
589 IFPACK_CHK_ERR(AminusAT.InsertGlobalValues(1,&GRID,1,&GCID,&plus_val));
592 if (AminusAT.SumIntoGlobalValues(1,&GCID,1,&GRID,&minus_val) != 0) {
593 IFPACK_CHK_ERR(AminusAT.InsertGlobalValues(1,&GCID,1,&GRID,&minus_val));
601 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 602 if(A.RowMatrixRowMap().GlobalIndicesLongLong()) {
603 for (
int i = 0 ; i < NumMyRows ; ++i) {
605 long long GRID = A.RowMatrixRowMap().GID64(i);
609 IFPACK_CHK_ERR(A.ExtractMyRowCopy(i,A.MaxNumEntries(),Nnz,
610 &colVal[0],&colInd[0]));
612 for (
int j = 0 ; j < Nnz ; ++j) {
614 long long GCID = A.RowMatrixColMap().GID64(colInd[j]);
617 double plus_val = colVal[j];
618 double minus_val = -colVal[j];
620 if (AplusAT.SumIntoGlobalValues(1,&GRID,1,&GCID,&plus_val) != 0) {
621 IFPACK_CHK_ERR(AplusAT.InsertGlobalValues(1,&GRID,1,&GCID,&plus_val));
624 if (AplusAT.SumIntoGlobalValues(1,&GCID,1,&GRID,&plus_val) != 0) {
625 IFPACK_CHK_ERR(AplusAT.InsertGlobalValues(1,&GCID,1,&GRID,&plus_val));
628 if (AminusAT.SumIntoGlobalValues(1,&GRID,1,&GCID,&plus_val) != 0) {
629 IFPACK_CHK_ERR(AminusAT.InsertGlobalValues(1,&GRID,1,&GCID,&plus_val));
632 if (AminusAT.SumIntoGlobalValues(1,&GCID,1,&GRID,&minus_val) != 0) {
633 IFPACK_CHK_ERR(AminusAT.InsertGlobalValues(1,&GCID,1,&GRID,&minus_val));
641 throw "Ifpack_Analyze: GlobalIndices type unknown";
643 AplusAT.FillComplete();
644 AminusAT.FillComplete();
649 NormOne = AplusAT.NormOne();
650 NormInf = AplusAT.NormInf();
651 NormF = Ifpack_FrobeniusNorm(AplusAT);
654 print<double>(
"A + A^T", NormOne, NormInf, NormF);
657 NormOne = AminusAT.NormOne();
658 NormInf = AminusAT.NormInf();
659 NormF = Ifpack_FrobeniusNorm(AminusAT);
662 print<double>(
"A - A^T", NormOne, NormInf, NormF);
668 print<const char*>(
"",
"min",
"avg",
"max",
false);
669 print<const char*>(
"",
"===",
"===",
"===",
false);
676 for (
int i = 0 ; i < NumMyRows ; ++i) {
679 IFPACK_CHK_ERR(A.ExtractMyRowCopy(i,A.MaxNumEntries(),Nnz,
680 &colVal[0],&colInd[0]));
682 for (
int j = 0 ; j < Nnz ; ++j) {
684 if (colVal[j] > MyMax) MyMax = colVal[j];
685 if (colVal[j] < MyMin) MyMin = colVal[j];
689 A.Comm().MaxAll(&MyMax, &GlobalMax, 1);
690 A.Comm().MinAll(&MyMin, &GlobalMin, 1);
691 A.Comm().SumAll(&MyAvg, &GlobalAvg, 1);
692 GlobalAvg /= A.NumGlobalNonzeros64();
696 print<double>(
" A(i,j)", GlobalMin, GlobalAvg, GlobalMax);
703 for (
int i = 0 ; i < NumMyRows ; ++i) {
706 IFPACK_CHK_ERR(A.ExtractMyRowCopy(i,A.MaxNumEntries(),Nnz,
707 &colVal[0],&colInd[0]));
709 for (
int j = 0 ; j < Nnz ; ++j) {
710 double v = colVal[j];
713 if (colVal[j] > MyMax) MyMax = v;
714 if (colVal[j] < MyMin) MyMin = v;
718 A.Comm().MaxAll(&MyMax, &GlobalMax, 1);
719 A.Comm().MinAll(&MyMin, &GlobalMin, 1);
720 A.Comm().SumAll(&MyAvg, &GlobalAvg, 1);
721 GlobalAvg /= A.NumGlobalNonzeros64();
724 print<double>(
"|A(i,j)|", GlobalMin, GlobalAvg, GlobalMax);
731 Diag.MinValue(&GlobalMin);
732 Diag.MaxValue(&GlobalMax);
733 Diag.MeanValue(&GlobalAvg);
737 print<double>(
" A(k,k)", GlobalMin, GlobalAvg, GlobalMax);
741 Diag.MinValue(&GlobalMin);
742 Diag.MaxValue(&GlobalMax);
743 Diag.MeanValue(&GlobalAvg);
745 print<double>(
"|A(k,k)|", GlobalMin, GlobalAvg, GlobalMax);
752 if (NumPDEEqns > 1 ) {
754 if (verbose) print();
756 for (
int ie = 0 ; ie < NumPDEEqns ; ie++) {
762 for (
int i = ie ; i < Diag.MyLength() ; i += NumPDEEqns) {
770 A.Comm().MinAll(&MyMin, &GlobalMin, 1);
771 A.Comm().MaxAll(&MyMax, &GlobalMax, 1);
772 A.Comm().SumAll(&MyAvg, &GlobalAvg, 1);
775 GlobalAvg /= (Diag.GlobalLength64() / NumPDEEqns);
779 sprintf(str,
" A(k,k), eq %d", ie);
780 print<double>(str, GlobalMin, GlobalAvg, GlobalMax);
789 RowSum.MinValue(&GlobalMin);
790 RowSum.MaxValue(&GlobalMax);
791 RowSum.MeanValue(&GlobalAvg);
795 print<double>(
" sum_j A(k,j)", GlobalMin, GlobalAvg, GlobalMax);
802 if (NumPDEEqns > 1 ) {
804 if (verbose) print();
806 for (
int ie = 0 ; ie < NumPDEEqns ; ie++) {
812 for (
int i = ie ; i < Diag.MyLength() ; i += NumPDEEqns) {
813 double d = RowSum[i];
820 A.Comm().MinAll(&MyMin, &GlobalMin, 1);
821 A.Comm().MaxAll(&MyMax, &GlobalMax, 1);
822 A.Comm().SumAll(&MyAvg, &GlobalAvg, 1);
825 GlobalAvg /= (Diag.GlobalLength64() / NumPDEEqns);
829 sprintf(str,
" sum_j A(k,j), eq %d", ie);
830 print<double>(str, GlobalMin, GlobalAvg, GlobalMax);
842 const bool abs,
const int steps)
847 bool verbose = (Diagonal.Comm().MyPID() == 0);
848 double min_val = DBL_MAX;
849 double max_val = -DBL_MAX;
851 for (
int i = 0 ; i < Diagonal.MyLength() ; ++i) {
852 double v = Diagonal[i];
864 cout <<
"Vector label = " << Diagonal.Label() << endl;
868 double delta = (max_val - min_val) / steps;
869 for (
int k = 0 ; k < steps ; ++k) {
871 double below = delta * k + min_val;
872 double above = below + delta;
873 int MyBelow = 0, GlobalBelow;
875 for (
int i = 0 ; i < Diagonal.MyLength() ; ++i) {
876 double v = Diagonal[i];
878 if (v >= below && v < above) MyBelow++;
881 Diagonal.Comm().SumAll(&MyBelow, &GlobalBelow, 1);
884 printf(
"Elements in [%+7e, %+7e) = %10d ( = %5.2f %%)\n",
885 below, above, GlobalBelow,
886 100.0 * GlobalBelow / Diagonal.GlobalLength64());
901 const bool abs,
const int steps)
906 bool verbose = (A.Comm().MyPID() == 0);
907 double min_val = DBL_MAX;
908 double max_val = -DBL_MAX;
910 std::vector<int> colInd(A.MaxNumEntries());
911 std::vector<double> colVal(A.MaxNumEntries());
913 for (
int i = 0 ; i < A.NumMyRows() ; ++i) {
916 IFPACK_CHK_ERR(A.ExtractMyRowCopy(i,A.MaxNumEntries(),Nnz,
917 &colVal[0],&colInd[0]));
919 for (
int j = 0 ; j < Nnz ; ++j) {
920 double v = colVal[j];
933 cout <<
"Label of matrix = " << A.Label() << endl;
937 double delta = (max_val - min_val) / steps;
938 for (
int k = 0 ; k < steps ; ++k) {
940 double below = delta * k + min_val;
941 double above = below + delta;
942 int MyBelow = 0, GlobalBelow;
944 for (
int i = 0 ; i < A.NumMyRows() ; ++i) {
947 IFPACK_CHK_ERR(A.ExtractMyRowCopy(i,A.MaxNumEntries(),Nnz,
948 &colVal[0],&colInd[0]));
950 for (
int j = 0 ; j < Nnz ; ++j) {
951 double v = colVal[j];
954 if (v >= below && v < above) MyBelow++;
958 A.Comm().SumAll(&MyBelow, &GlobalBelow, 1);
960 printf(
"Elements in [%+7e, %+7e) = %10d ( = %5.2f %%)\n",
961 below, above, GlobalBelow,
962 100.0 * GlobalBelow / A.NumGlobalNonzeros64());
975 int Ifpack_PrintSparsity(
const Epetra_RowMatrix& A,
const char* InputFileName,
976 const int NumPDEEqns)
980 long long m,nc,nr,maxdim;
981 double lrmrgn,botmrgn,xtit,ytit,ytitof,fnstit,siz = 0.0;
982 double xl,xr, yb,yt, scfct,u2dot,frlw,delt,paperx;
993 long long NumGlobalRows;
994 long long NumGlobalCols;
1000 const Epetra_Comm& Comm = A.Comm();
1004 if (strlen(A.Label()) != 0)
1005 strcpy(title, A.Label());
1007 sprintf(title,
"%s",
"matrix");
1009 if (InputFileName == 0)
1010 sprintf(FileName,
"%s.ps", title);
1012 strcpy(FileName, InputFileName);
1014 MyPID = Comm.MyPID();
1015 NumProc = Comm.NumProc();
1017 NumMyRows = A.NumMyRows();
1020 NumGlobalRows = A.NumGlobalRows64();
1021 NumGlobalCols = A.NumGlobalCols64();
1023 if (NumGlobalRows != NumGlobalCols)
1027 maxdim = (NumGlobalRows>NumGlobalCols)?NumGlobalRows:NumGlobalCols;
1028 maxdim /= NumPDEEqns;
1031 nr = NumGlobalRows / NumPDEEqns + 1;
1032 nc = NumGlobalCols / NumPDEEqns + 1;
1047 lrmrgn = (paperx-siz)/2.0;
1053 scfct = siz*u2dot/m;
1063 ltit = strlen(title);
1069 ytit = botmrgn+siz*nr/m + ytitof;
1071 xl = lrmrgn*u2dot - scfct*frlw/2;
1072 xr = (lrmrgn+siz)*u2dot + scfct*frlw/2;
1073 yb = botmrgn*u2dot - scfct*frlw/2;
1074 yt = (botmrgn+siz*nr/m)*u2dot + scfct*frlw/2;
1076 yt = yt + (ytitof+fnstit*0.70)*u2dot;
1086 if ((ptitle == 0) && (ltit == 0)) {
1087 ytit = botmrgn + fnstit*0.3;
1088 botmrgn = botmrgn + ytitof + fnstit*0.7;
1095 fp = fopen(FileName,
"w");
1097 fprintf(fp,
"%s",
"%%!PS-Adobe-2.0\n");
1098 fprintf(fp,
"%s",
"%%Creator: IFPACK\n");
1099 fprintf(fp,
"%%%%BoundingBox: %f %f %f %f\n",
1101 fprintf(fp,
"%s",
"%%EndComments\n");
1102 fprintf(fp,
"%s",
"/cm {72 mul 2.54 div} def\n");
1103 fprintf(fp,
"%s",
"/mc {72 div 2.54 mul} def\n");
1104 fprintf(fp,
"%s",
"/pnum { 72 div 2.54 mul 20 std::string ");
1105 fprintf(fp,
"%s",
"cvs print ( ) print} def\n");
1106 fprintf(fp,
"%s",
"/Cshow {dup stringwidth pop -2 div 0 rmoveto show} def\n");
1110 fprintf(fp,
"%s",
"gsave\n");
1112 fprintf(fp,
"/Helvetica findfont %e cm scalefont setfont\n",
1114 fprintf(fp,
"%f cm %f cm moveto\n",
1116 fprintf(fp,
"(%s) Cshow\n", title);
1117 fprintf(fp,
"%f cm %f cm translate\n",
1120 fprintf(fp,
"%f cm %d div dup scale \n",
1124 fprintf(fp,
"%f setlinewidth\n",
1126 fprintf(fp,
"%s",
"newpath\n");
1127 fprintf(fp,
"%s",
"0 0 moveto ");
1129 printf(
"------------------- %d\n", (
int) m);
1130 fprintf(fp,
"%d %d lineto\n",
1132 fprintf(fp,
"%d %d lineto\n",
1134 fprintf(fp,
"%d %d lineto\n",
1138 fprintf(fp,
"%d %d lineto\n",
1140 fprintf(fp,
"%d %d lineto\n",
1141 (
int) nc, (
int) nr);
1142 fprintf(fp,
"%d %d lineto\n",
1145 fprintf(fp,
"%s",
"closepath stroke\n");
1149 fprintf(fp,
"%s",
"1 1 translate\n");
1150 fprintf(fp,
"%s",
"0.8 setlinewidth\n");
1151 fprintf(fp,
"%s",
"/p {moveto 0 -.40 rmoveto \n");
1152 fprintf(fp,
"%s",
" 0 .80 rlineto stroke} def\n");
1157 int MaxEntries = A.MaxNumEntries();
1158 std::vector<int> Indices(MaxEntries);
1159 std::vector<double> Values(MaxEntries);
1161 for (
int pid = 0 ; pid < NumProc ; ++pid) {
1165 fp = fopen(FileName,
"a");
1166 TEUCHOS_ASSERT(fp != NULL);
1168 for (
int i = 0 ; i < NumMyRows ; ++i) {
1170 if (i % NumPDEEqns)
continue;
1173 A.ExtractMyRowCopy(i,MaxEntries,Nnz,&Values[0],&Indices[0]);
1175 long long grow = A.RowMatrixRowMap().GID64(i);
1177 for (
int j = 0 ; j < Nnz ; ++j) {
1178 int col = Indices[j];
1179 if (col % NumPDEEqns == 0) {
1180 long long gcol = A.RowMatrixColMap().GID64(Indices[j]);
1183 fprintf(fp,
"%lld %lld p\n",
1184 gcol, NumGlobalRows - grow - 1);
1189 fprintf(fp,
"%s",
"%end of data for this process\n");
1191 if( pid == NumProc - 1 )
1192 fprintf(fp,
"%s",
"showpage\n");
int Ifpack_Analyze(const Epetra_RowMatrix &A, const bool Cheap=false, const int NumPDEEqns=1)
Analyzes the basic properties of the input matrix A; see Usage of Ifpack_Analyze()..
int Ifpack_AnalyzeVectorElements(const Epetra_Vector &Diagonal, const bool abs=false, const int steps=10)
Analyzes the distribution of values of the input vector Diagonal.
int Ifpack_AnalyzeMatrixElements(const Epetra_RowMatrix &A, const bool abs=false, const int steps=10)
Analyzes the distribution of values of the input matrix A.
void Ifpack_PrintLine()
Prints a line of ‘=’ on cout.
Epetra_CrsMatrix * Ifpack_CreateOverlappingCrsMatrix(const Epetra_RowMatrix *Matrix, const int OverlappingLevel)
Creates an overlapping Epetra_CrsMatrix. Returns 0 if OverlappingLevel is 0.
void Ifpack_BreakForDebugger(Epetra_Comm &Comm)
Stops the execution of code, so that a debugger can be attached.
int Ifpack_PrintResidual(char *Label, const Epetra_RowMatrix &A, const Epetra_MultiVector &X, const Epetra_MultiVector &Y)
Prints on cout the true residual.
std::string Ifpack_toString(const int &x)
Converts an integer to std::string.