IFPACK Development
Loading...
Searching...
No Matches
Ifpack_Utils.cpp
1/*@HEADER
2// ***********************************************************************
3//
4// Ifpack: Object-Oriented Algebraic Preconditioner Package
5// Copyright (2002) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38//
39// ***********************************************************************
40//@HEADER
41*/
42
43#include "Ifpack_ConfigDefs.h"
44#include "Ifpack_Preconditioner.h"
45#include "Ifpack_Utils.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"
54
55void Ifpack_PrintLine()
56{
57 std::cout << "================================================================================" << std::endl;
58}
59
60//============================================================================
61void Ifpack_BreakForDebugger(Epetra_Comm& Comm)
62{
63 using std::cout;
64 using std::endl;
65
66 char hostname[80];
67 char buf[160];
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());
73#elif defined(_WIN32)
74 sprintf(buf,"Windows compiler, unknown hostname and PID!");
75#else
76 gethostname(hostname, sizeof(hostname));
77 sprintf(buf, "Host: %s\tComm.MyPID(): %d\tPID: %d",
78 hostname, Comm.MyPID(), getpid());
79#endif
80 printf("%s\n",buf);
81 fflush(stdout);
82#if !( defined(_WIN32) )
83 sleep(1);
84#endif
85 }
86 }
87 if(Comm.MyPID() == 0) {
88 printf("\n");
89 printf("** Pausing to attach debugger...\n");
90 printf("** You may now attach debugger to the processes listed above.\n");
91 printf( "**\n");
92 printf( "** Enter a character to continue > "); fflush(stdout);
93 char go;
94 TEUCHOS_ASSERT(scanf("%c",&go) != EOF);
95 }
96
97 Comm.Barrier();
98
99}
100
101//============================================================================
102Epetra_CrsMatrix* Ifpack_CreateOverlappingCrsMatrix(const Epetra_RowMatrix* Matrix,
103 const int OverlappingLevel)
104{
105
106 if (OverlappingLevel == 0)
107 return(0); // All done
108 if (Matrix->Comm().NumProc() == 1)
109 return(0); // All done
110
111 Epetra_CrsMatrix* OverlappingMatrix;
112 OverlappingMatrix = 0;
113 Epetra_Map* OverlappingMap;
114 OverlappingMap = (Epetra_Map*)&(Matrix->RowMatrixRowMap());
115
116 const Epetra_RowMatrix* OldMatrix;
117 const Epetra_Map* DomainMap = &(Matrix->OperatorDomainMap());
118 const Epetra_Map* RangeMap = &(Matrix->OperatorRangeMap());
119
120 for (int level = 1; level <= OverlappingLevel ; ++level) {
121
122 if (OverlappingMatrix)
123 OldMatrix = OverlappingMatrix;
124 else
125 OldMatrix = Matrix;
126
127 Epetra_Import* OverlappingImporter;
128 OverlappingImporter = (Epetra_Import*)OldMatrix->RowMatrixImporter();
129 int NumMyElements = OverlappingImporter->TargetMap().NumMyElements();
130
131 // need to build an Epetra_Map in this way because Epetra_CrsMatrix
132 // requires Epetra_Map and not Epetra_BlockMap
133
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,
138 0, Matrix->Comm());
139 }
140 else
141#endif
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,
146 0, Matrix->Comm());
147 }
148 else
149#endif
150 throw "Ifpack_CreateOverlappingCrsMatrix: GlobalIndices type unknown";
151
152 if (level < OverlappingLevel)
153 OverlappingMatrix = new Epetra_CrsMatrix(Copy, *OverlappingMap, 0);
154 else
155 // On last iteration, we want to filter out all columns except
156 // those that correspond
157 // to rows in the graph. This assures that our matrix is square
158 OverlappingMatrix = new Epetra_CrsMatrix(Copy, *OverlappingMap,
159 *OverlappingMap, 0);
160
161 OverlappingMatrix->Import(*OldMatrix, *OverlappingImporter, Insert);
162 if (level < OverlappingLevel) {
163 OverlappingMatrix->FillComplete(*DomainMap, *RangeMap);
164 }
165 else {
166 OverlappingMatrix->FillComplete(*DomainMap, *RangeMap);
167 }
168
169 delete OverlappingMap;
170
171 if (level > 1) {
172 delete OldMatrix;
173 }
174 OverlappingMatrix->FillComplete();
175
176 }
177
178 return(OverlappingMatrix);
179}
180
181//============================================================================
182Epetra_CrsGraph* Ifpack_CreateOverlappingCrsMatrix(const Epetra_CrsGraph* Graph,
183 const int OverlappingLevel)
184{
185
186 if (OverlappingLevel == 0)
187 return(0); // All done
188 if (Graph->Comm().NumProc() == 1)
189 return(0); // All done
190
191 Epetra_CrsGraph* OverlappingGraph;
192 Epetra_BlockMap* OverlappingMap;
193 OverlappingGraph = const_cast<Epetra_CrsGraph*>(Graph);
194 OverlappingMap = const_cast<Epetra_BlockMap*>(&(Graph->RowMap()));
195
196 Epetra_CrsGraph* OldGraph;
197 Epetra_BlockMap* OldMap;
198 const Epetra_BlockMap* DomainMap = &(Graph->DomainMap());
199 const Epetra_BlockMap* RangeMap = &(Graph->RangeMap());
200
201 for (int level = 1; level <= OverlappingLevel ; ++level) {
202
203 OldGraph = OverlappingGraph;
204 OldMap = OverlappingMap;
205
206 Epetra_Import* OverlappingImporter;
207 OverlappingImporter = const_cast<Epetra_Import*>(OldGraph->Importer());
208 OverlappingMap = new Epetra_BlockMap(OverlappingImporter->TargetMap());
209
210 if (level < OverlappingLevel)
211 OverlappingGraph = new Epetra_CrsGraph(Copy, *OverlappingMap, 0);
212 else
213 // On last iteration, we want to filter out all columns except
214 // those that correspond
215 // to rows in the graph. This assures that our matrix is square
216 OverlappingGraph = new Epetra_CrsGraph(Copy, *OverlappingMap,
217 *OverlappingMap, 0);
218
219 OverlappingGraph->Import(*OldGraph, *OverlappingImporter, Insert);
220 if (level < OverlappingLevel)
221 OverlappingGraph->FillComplete(*DomainMap, *RangeMap);
222 else {
223 // Copy last OverlapImporter because we will use it later
224 OverlappingImporter = new Epetra_Import(*OverlappingMap, *DomainMap);
225 OverlappingGraph->FillComplete(*DomainMap, *RangeMap);
226 }
227
228 if (level > 1) {
229 delete OldGraph;
230 delete OldMap;
231 }
232
233 delete OverlappingMap;
234 OverlappingGraph->FillComplete();
235
236 }
237
238 return(OverlappingGraph);
239}
240
241//============================================================================
242std::string Ifpack_toString(const int& x)
243{
244 char s[100];
245 sprintf(s, "%d", x);
246 return std::string(s);
247}
248
249//============================================================================
250std::string Ifpack_toString(const double& x)
251{
252 char s[100];
253 sprintf(s, "%g", x);
254 return std::string(s);
255}
256
257//============================================================================
258int Ifpack_PrintResidual(char* Label, const Epetra_RowMatrix& A,
259 const Epetra_MultiVector& X, const Epetra_MultiVector&Y)
260{
261 using std::cout;
262 using std::endl;
263
264 if (X.Comm().MyPID() == 0) {
265 cout << "***** " << Label << endl;
266 }
267 Ifpack_PrintResidual(0,A,X,Y);
268
269 return(0);
270}
271
272//============================================================================
273int Ifpack_PrintResidual(const int iter, const Epetra_RowMatrix& A,
274 const Epetra_MultiVector& X, const Epetra_MultiVector&Y)
275{
276 using std::cout;
277 using std::endl;
278
279 Epetra_MultiVector RHS(X);
280 std::vector<double> Norm2;
281 Norm2.resize(X.NumVectors());
282
283 IFPACK_CHK_ERR(A.Multiply(false,X,RHS));
284 RHS.Update(1.0, Y, -1.0);
285
286 RHS.Norm2(&Norm2[0]);
287
288 if (X.Comm().MyPID() == 0) {
289 cout << "***** iter: " << iter << ": ||Ax - b||_2 = "
290 << Norm2[0] << endl;
291 }
292
293 return(0);
294}
295
296//============================================================================
297// very simple debugging function that prints on screen the structure
298// of the input matrix.
299void Ifpack_PrintSparsity_Simple(const Epetra_RowMatrix& A)
300{
301 using std::cout;
302 using std::endl;
303
304 int MaxEntries = A.MaxNumEntries();
305 std::vector<int> Indices(MaxEntries);
306 std::vector<double> Values(MaxEntries);
307 std::vector<bool> FullRow(A.NumMyRows());
308
309 cout << "+-";
310 for (int j = 0 ; j < A.NumMyRows() ; ++j)
311 cout << '-';
312 cout << "-+" << endl;
313
314 for (int i = 0 ; i < A.NumMyRows() ; ++i) {
315
316 int Length;
317 A.ExtractMyRowCopy(i,MaxEntries,Length,
318 &Values[0], &Indices[0]);
319
320 for (int j = 0 ; j < A.NumMyRows() ; ++j)
321 FullRow[j] = false;
322
323 for (int j = 0 ; j < Length ; ++j) {
324 FullRow[Indices[j]] = true;
325 }
326
327 cout << "| ";
328 for (int j = 0 ; j < A.NumMyRows() ; ++j) {
329 if (FullRow[j])
330 cout << '*';
331 else
332 cout << ' ';
333 }
334 cout << " |" << endl;
335 }
336
337 cout << "+-";
338 for (int j = 0 ; j < A.NumMyRows() ; ++j)
339 cout << '-';
340 cout << "-+" << endl << endl;
341
342}
343
344//============================================================================
345
346double Ifpack_FrobeniusNorm(const Epetra_RowMatrix& A)
347{
348 double MyNorm = 0.0, GlobalNorm;
349
350 std::vector<int> colInd(A.MaxNumEntries());
351 std::vector<double> colVal(A.MaxNumEntries());
352
353 for (int i = 0 ; i < A.NumMyRows() ; ++i) {
354
355 int Nnz;
356 IFPACK_CHK_ERR(A.ExtractMyRowCopy(i,A.MaxNumEntries(),Nnz,
357 &colVal[0],&colInd[0]));
358
359 for (int j = 0 ; j < Nnz ; ++j) {
360 MyNorm += colVal[j] * colVal[j];
361 }
362 }
363
364 A.Comm().SumAll(&MyNorm,&GlobalNorm,1);
365
366 return(sqrt(GlobalNorm));
367}
368
369static void print()
370{
371 printf("\n");
372}
373
374#include <iomanip>
375template<class T>
376static void print(const char str[], T val)
377{
378 std::cout.width(30); std::cout.setf(std::ios::left);
379 std::cout << str;
380 std::cout << " = " << val << std::endl;
381}
382
383template<class T>
384static void print(const char str[], T val, double percentage)
385{
386 std::cout.width(30); std::cout.setf(std::ios::left);
387 std::cout << str;
388 std::cout << " = ";
389 std::cout.width(20); std::cout.setf(std::ios::left);
390 std::cout << val;
391 std::cout << " ( " << percentage << " %)" << std::endl;
392}
393template<class T>
394static void print(const char str[], T one, T two, T three, bool equal = true)
395{
396 using std::endl;
397
398 std::cout.width(30); std::cout.setf(std::ios::left);
399 std::cout << str;
400 if (equal)
401 std::cout << " = ";
402 else
403 std::cout << " ";
404 std::cout.width(15); std::cout.setf(std::ios::left);
405 std::cout << one;
406 std::cout.width(15); std::cout.setf(std::ios::left);
407 std::cout << two;
408 std::cout.width(15); std::cout.setf(std::ios::left);
409 std::cout << three;
410 std::cout << endl;
411}
412
413//============================================================================
414#include "limits.h"
415#include "float.h"
416#include "Epetra_FECrsMatrix.h"
417
418int Ifpack_Analyze(const Epetra_RowMatrix& A, const bool Cheap,
419 const int NumPDEEqns)
420{
421
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;
433
434 bool verbose = (A.Comm().MyPID() == 0);
435
436 GlobalStorage = sizeof(int*) * NumGlobalRows +
437 sizeof(int) * A.NumGlobalNonzeros64() +
438 sizeof(double) * A.NumGlobalNonzeros64();
439
440 if (verbose) {
441 print();
442 Ifpack_PrintLine();
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);
449 }
450
451 long long NumMyActualNonzeros = 0, NumGlobalActualNonzeros;
452 long long NumMyEmptyRows = 0, NumGlobalEmptyRows;
453 long long NumMyDirichletRows = 0, NumGlobalDirichletRows;
454
455 std::vector<int> colInd(A.MaxNumEntries());
456 std::vector<double> colVal(A.MaxNumEntries());
457
459 Epetra_Vector RowSum(A.RowMatrixRowMap());
460 Diag.PutScalar(0.0);
461 RowSum.PutScalar(0.0);
462
463 for (int i = 0 ; i < NumMyRows ; ++i) {
464
465 long long GRID = A.RowMatrixRowMap().GID64(i);
466 int Nnz;
467 IFPACK_CHK_ERR(A.ExtractMyRowCopy(i,A.MaxNumEntries(),Nnz,
468 &colVal[0],&colInd[0]));
469
470 if (Nnz == 0)
471 NumMyEmptyRows++;
472
473 if (Nnz == 1)
474 NumMyDirichletRows++;
475
476 for (int j = 0 ; j < Nnz ; ++j) {
477
478 double v = colVal[j];
479 if (v < 0) v = -v;
480 if (colVal[j] != 0.0)
481 NumMyActualNonzeros++;
482
483 long long GCID = A.RowMatrixColMap().GID64(colInd[j]);
484
485 if (GCID != GRID)
486 RowSum[i] += v;
487 else
488 Diag[i] = v;
489
490 if (GCID < GRID)
491 MyLowerNonzeros++;
492 else if (GCID > GRID)
493 MyUpperNonzeros++;
494 long long b = GCID - GRID;
495 if (b < 0) b = -b;
496 if (b > MyBandwidth)
497 MyBandwidth = b;
498 }
499
500 if (Diag[i] > RowSum[i])
501 MyDiagonallyDominant++;
502
503 if (Diag[i] >= RowSum[i])
504 MyWeaklyDiagonallyDominant++;
505
506 RowSum[i] += Diag[i];
507 }
508
509 // ======================== //
510 // summing up global values //
511 // ======================== //
512
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);
523
524 double NormOne = A.NormOne();
525 double NormInf = A.NormInf();
526 double NormF = Ifpack_FrobeniusNorm(A);
527
528 if (verbose) {
529 print();
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);
533 print();
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);
543 print();
544 print<long long>("Maximum bandwidth", GlobalBandwidth);
545
546 print();
547 print("", "one-norm", "inf-norm", "Frobenius", false);
548 print("", "========", "========", "=========", false);
549 print();
550
551 print<double>("A", NormOne, NormInf, NormF);
552 }
553
554 if (Cheap == false) {
555
556 // create A + A^T and A - A^T
557
558 Epetra_FECrsMatrix AplusAT(Copy, A.RowMatrixRowMap(), 0);
559 Epetra_FECrsMatrix AminusAT(Copy, A.RowMatrixRowMap(), 0);
560
561#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
563 for (int i = 0 ; i < NumMyRows ; ++i) {
564
565 int GRID = A.RowMatrixRowMap().GID(i);
566 assert (GRID != -1);
567
568 int Nnz;
569 IFPACK_CHK_ERR(A.ExtractMyRowCopy(i,A.MaxNumEntries(),Nnz,
570 &colVal[0],&colInd[0]));
571
572 for (int j = 0 ; j < Nnz ; ++j) {
573
574 int GCID = A.RowMatrixColMap().GID(colInd[j]);
575 assert (GCID != -1);
576
577 double plus_val = colVal[j];
578 double minus_val = -colVal[j];
579
580 if (AplusAT.SumIntoGlobalValues(1,&GRID,1,&GCID,&plus_val) != 0) {
581 IFPACK_CHK_ERR(AplusAT.InsertGlobalValues(1,&GRID,1,&GCID,&plus_val));
582 }
583
584 if (AplusAT.SumIntoGlobalValues(1,&GCID,1,&GRID,&plus_val) != 0) {
585 IFPACK_CHK_ERR(AplusAT.InsertGlobalValues(1,&GCID,1,&GRID,&plus_val));
586 }
587
588 if (AminusAT.SumIntoGlobalValues(1,&GRID,1,&GCID,&plus_val) != 0) {
589 IFPACK_CHK_ERR(AminusAT.InsertGlobalValues(1,&GRID,1,&GCID,&plus_val));
590 }
591
592 if (AminusAT.SumIntoGlobalValues(1,&GCID,1,&GRID,&minus_val) != 0) {
593 IFPACK_CHK_ERR(AminusAT.InsertGlobalValues(1,&GCID,1,&GRID,&minus_val));
594 }
595
596 }
597 }
598 }
599 else
600#endif
601#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
603 for (int i = 0 ; i < NumMyRows ; ++i) {
604
605 long long GRID = A.RowMatrixRowMap().GID64(i);
606 assert (GRID != -1);
607
608 int Nnz;
609 IFPACK_CHK_ERR(A.ExtractMyRowCopy(i,A.MaxNumEntries(),Nnz,
610 &colVal[0],&colInd[0]));
611
612 for (int j = 0 ; j < Nnz ; ++j) {
613
614 long long GCID = A.RowMatrixColMap().GID64(colInd[j]);
615 assert (GCID != -1);
616
617 double plus_val = colVal[j];
618 double minus_val = -colVal[j];
619
620 if (AplusAT.SumIntoGlobalValues(1,&GRID,1,&GCID,&plus_val) != 0) {
621 IFPACK_CHK_ERR(AplusAT.InsertGlobalValues(1,&GRID,1,&GCID,&plus_val));
622 }
623
624 if (AplusAT.SumIntoGlobalValues(1,&GCID,1,&GRID,&plus_val) != 0) {
625 IFPACK_CHK_ERR(AplusAT.InsertGlobalValues(1,&GCID,1,&GRID,&plus_val));
626 }
627
628 if (AminusAT.SumIntoGlobalValues(1,&GRID,1,&GCID,&plus_val) != 0) {
629 IFPACK_CHK_ERR(AminusAT.InsertGlobalValues(1,&GRID,1,&GCID,&plus_val));
630 }
631
632 if (AminusAT.SumIntoGlobalValues(1,&GCID,1,&GRID,&minus_val) != 0) {
633 IFPACK_CHK_ERR(AminusAT.InsertGlobalValues(1,&GCID,1,&GRID,&minus_val));
634 }
635
636 }
637 }
638 }
639 else
640#endif
641 throw "Ifpack_Analyze: GlobalIndices type unknown";
642
643 AplusAT.FillComplete();
644 AminusAT.FillComplete();
645
646 AplusAT.Scale(0.5);
647 AminusAT.Scale(0.5);
648
649 NormOne = AplusAT.NormOne();
650 NormInf = AplusAT.NormInf();
651 NormF = Ifpack_FrobeniusNorm(AplusAT);
652
653 if (verbose) {
654 print<double>("A + A^T", NormOne, NormInf, NormF);
655 }
656
657 NormOne = AminusAT.NormOne();
658 NormInf = AminusAT.NormInf();
659 NormF = Ifpack_FrobeniusNorm(AminusAT);
660
661 if (verbose) {
662 print<double>("A - A^T", NormOne, NormInf, NormF);
663 }
664 }
665
666 if (verbose) {
667 print();
668 print<const char*>("", "min", "avg", "max", false);
669 print<const char*>("", "===", "===", "===", false);
670 }
671
672 MyMax = -DBL_MAX;
673 MyMin = DBL_MAX;
674 MyAvg = 0.0;
675
676 for (int i = 0 ; i < NumMyRows ; ++i) {
677
678 int Nnz;
679 IFPACK_CHK_ERR(A.ExtractMyRowCopy(i,A.MaxNumEntries(),Nnz,
680 &colVal[0],&colInd[0]));
681
682 for (int j = 0 ; j < Nnz ; ++j) {
683 MyAvg += colVal[j];
684 if (colVal[j] > MyMax) MyMax = colVal[j];
685 if (colVal[j] < MyMin) MyMin = colVal[j];
686 }
687 }
688
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();
693
694 if (verbose) {
695 print();
696 print<double>(" A(i,j)", GlobalMin, GlobalAvg, GlobalMax);
697 }
698
699 MyMax = 0.0;
700 MyMin = DBL_MAX;
701 MyAvg = 0.0;
702
703 for (int i = 0 ; i < NumMyRows ; ++i) {
704
705 int Nnz;
706 IFPACK_CHK_ERR(A.ExtractMyRowCopy(i,A.MaxNumEntries(),Nnz,
707 &colVal[0],&colInd[0]));
708
709 for (int j = 0 ; j < Nnz ; ++j) {
710 double v = colVal[j];
711 if (v < 0) v = -v;
712 MyAvg += v;
713 if (colVal[j] > MyMax) MyMax = v;
714 if (colVal[j] < MyMin) MyMin = v;
715 }
716 }
717
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();
722
723 if (verbose) {
724 print<double>("|A(i,j)|", GlobalMin, GlobalAvg, GlobalMax);
725 }
726
727 // ================= //
728 // diagonal elements //
729 // ================= //
730
731 Diag.MinValue(&GlobalMin);
732 Diag.MaxValue(&GlobalMax);
733 Diag.MeanValue(&GlobalAvg);
734
735 if (verbose) {
736 print();
737 print<double>(" A(k,k)", GlobalMin, GlobalAvg, GlobalMax);
738 }
739
740 Diag.Abs(Diag);
741 Diag.MinValue(&GlobalMin);
742 Diag.MaxValue(&GlobalMax);
743 Diag.MeanValue(&GlobalAvg);
744 if (verbose) {
745 print<double>("|A(k,k)|", GlobalMin, GlobalAvg, GlobalMax);
746 }
747
748 // ============================================== //
749 // cycle over all equations for diagonal elements //
750 // ============================================== //
751
752 if (NumPDEEqns > 1 ) {
753
754 if (verbose) print();
755
756 for (int ie = 0 ; ie < NumPDEEqns ; ie++) {
757
758 MyMin = DBL_MAX;
759 MyMax = -DBL_MAX;
760 MyAvg = 0.0;
761
762 for (int i = ie ; i < Diag.MyLength() ; i += NumPDEEqns) {
763 double d = Diag[i];
764 MyAvg += d;
765 if (d < MyMin)
766 MyMin = d;
767 if (d > MyMax)
768 MyMax = d;
769 }
770 A.Comm().MinAll(&MyMin, &GlobalMin, 1);
771 A.Comm().MaxAll(&MyMax, &GlobalMax, 1);
772 A.Comm().SumAll(&MyAvg, &GlobalAvg, 1);
773 // does not really work fine if the number of global
774 // elements is not a multiple of NumPDEEqns
775 GlobalAvg /= (Diag.GlobalLength64() / NumPDEEqns);
776
777 if (verbose) {
778 char str[80];
779 sprintf(str, " A(k,k), eq %d", ie);
780 print<double>(str, GlobalMin, GlobalAvg, GlobalMax);
781 }
782 }
783 }
784
785 // ======== //
786 // row sums //
787 // ======== //
788
789 RowSum.MinValue(&GlobalMin);
790 RowSum.MaxValue(&GlobalMax);
791 RowSum.MeanValue(&GlobalAvg);
792
793 if (verbose) {
794 print();
795 print<double>(" sum_j A(k,j)", GlobalMin, GlobalAvg, GlobalMax);
796 }
797
798 // ===================================== //
799 // cycle over all equations for row sums //
800 // ===================================== //
801
802 if (NumPDEEqns > 1 ) {
803
804 if (verbose) print();
805
806 for (int ie = 0 ; ie < NumPDEEqns ; ie++) {
807
808 MyMin = DBL_MAX;
809 MyMax = -DBL_MAX;
810 MyAvg = 0.0;
811
812 for (int i = ie ; i < Diag.MyLength() ; i += NumPDEEqns) {
813 double d = RowSum[i];
814 MyAvg += d;
815 if (d < MyMin)
816 MyMin = d;
817 if (d > MyMax)
818 MyMax = d;
819 }
820 A.Comm().MinAll(&MyMin, &GlobalMin, 1);
821 A.Comm().MaxAll(&MyMax, &GlobalMax, 1);
822 A.Comm().SumAll(&MyAvg, &GlobalAvg, 1);
823 // does not really work fine if the number of global
824 // elements is not a multiple of NumPDEEqns
825 GlobalAvg /= (Diag.GlobalLength64() / NumPDEEqns);
826
827 if (verbose) {
828 char str[80];
829 sprintf(str, " sum_j A(k,j), eq %d", ie);
830 print<double>(str, GlobalMin, GlobalAvg, GlobalMax);
831 }
832 }
833 }
834
835 if (verbose)
836 Ifpack_PrintLine();
837
838 return(0);
839}
840
841int Ifpack_AnalyzeVectorElements(const Epetra_Vector& Diagonal,
842 const bool abs, const int steps)
843{
844 using std::cout;
845 using std::endl;
846
847 bool verbose = (Diagonal.Comm().MyPID() == 0);
848 double min_val = DBL_MAX;
849 double max_val = -DBL_MAX;
850
851 for (int i = 0 ; i < Diagonal.MyLength() ; ++i) {
852 double v = Diagonal[i];
853 if (abs)
854 if (v < 0) v = -v;
855 if (v > max_val)
856 max_val = v;
857 if (v < min_val)
858 min_val = v;
859 }
860
861 if (verbose) {
862 cout << endl;
863 Ifpack_PrintLine();
864 cout << "Vector label = " << Diagonal.Label() << endl;
865 cout << endl;
866 }
867
868 double delta = (max_val - min_val) / steps;
869 for (int k = 0 ; k < steps ; ++k) {
870
871 double below = delta * k + min_val;
872 double above = below + delta;
873 int MyBelow = 0, GlobalBelow;
874
875 for (int i = 0 ; i < Diagonal.MyLength() ; ++i) {
876 double v = Diagonal[i];
877 if (v < 0) v = -v;
878 if (v >= below && v < above) MyBelow++;
879 }
880
881 Diagonal.Comm().SumAll(&MyBelow, &GlobalBelow, 1);
882
883 if (verbose) {
884 printf("Elements in [%+7e, %+7e) = %10d ( = %5.2f %%)\n",
885 below, above, GlobalBelow,
886 100.0 * GlobalBelow / Diagonal.GlobalLength64());
887 }
888 }
889
890 if (verbose) {
891 Ifpack_PrintLine();
892 cout << endl;
893 }
894
895 return(0);
896}
897
898// ======================================================================
899
900int Ifpack_AnalyzeMatrixElements(const Epetra_RowMatrix& A,
901 const bool abs, const int steps)
902{
903 using std::cout;
904 using std::endl;
905
906 bool verbose = (A.Comm().MyPID() == 0);
907 double min_val = DBL_MAX;
908 double max_val = -DBL_MAX;
909
910 std::vector<int> colInd(A.MaxNumEntries());
911 std::vector<double> colVal(A.MaxNumEntries());
912
913 for (int i = 0 ; i < A.NumMyRows() ; ++i) {
914
915 int Nnz;
916 IFPACK_CHK_ERR(A.ExtractMyRowCopy(i,A.MaxNumEntries(),Nnz,
917 &colVal[0],&colInd[0]));
918
919 for (int j = 0 ; j < Nnz ; ++j) {
920 double v = colVal[j];
921 if (abs)
922 if (v < 0) v = -v;
923 if (v < min_val)
924 min_val = v;
925 if (v > max_val)
926 max_val = v;
927 }
928 }
929
930 if (verbose) {
931 cout << endl;
932 Ifpack_PrintLine();
933 cout << "Label of matrix = " << A.Label() << endl;
934 cout << endl;
935 }
936
937 double delta = (max_val - min_val) / steps;
938 for (int k = 0 ; k < steps ; ++k) {
939
940 double below = delta * k + min_val;
941 double above = below + delta;
942 int MyBelow = 0, GlobalBelow;
943
944 for (int i = 0 ; i < A.NumMyRows() ; ++i) {
945
946 int Nnz;
947 IFPACK_CHK_ERR(A.ExtractMyRowCopy(i,A.MaxNumEntries(),Nnz,
948 &colVal[0],&colInd[0]));
949
950 for (int j = 0 ; j < Nnz ; ++j) {
951 double v = colVal[j];
952 if (abs)
953 if (v < 0) v = -v;
954 if (v >= below && v < above) MyBelow++;
955 }
956
957 }
958 A.Comm().SumAll(&MyBelow, &GlobalBelow, 1);
959 if (verbose) {
960 printf("Elements in [%+7e, %+7e) = %10d ( = %5.2f %%)\n",
961 below, above, GlobalBelow,
962 100.0 * GlobalBelow / A.NumGlobalNonzeros64());
963 }
964 }
965
966 if (verbose) {
967 Ifpack_PrintLine();
968 cout << endl;
969 }
970
971 return(0);
972}
973
974// ======================================================================
975int Ifpack_PrintSparsity(const Epetra_RowMatrix& A, const char* InputFileName,
976 const int NumPDEEqns)
977{
978
979 int ltit;
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;
983 bool square = false;
984 /*change square to .true. if you prefer a square frame around
985 a rectangular matrix */
986 double conv = 2.54;
987 char munt = 'E'; /* put 'E' for centimeters, 'U' for inches */
988 int ptitle = 0; /* position of the title, 0 under the drawing,
989 else above */
990 FILE* fp = NULL;
991 int NumMyRows;
992 //int NumMyCols;
993 long long NumGlobalRows;
994 long long NumGlobalCols;
995 int MyPID;
996 int NumProc;
997 char FileName[1024];
998 char title[1020];
999
1000 const Epetra_Comm& Comm = A.Comm();
1001
1002 /* --------------------- execution begins ---------------------- */
1003
1004 if (strlen(A.Label()) != 0)
1005 strcpy(title, A.Label());
1006 else
1007 sprintf(title, "%s", "matrix");
1008
1009 if (InputFileName == 0)
1010 sprintf(FileName, "%s.ps", title);
1011 else
1012 strcpy(FileName, InputFileName);
1013
1014 MyPID = Comm.MyPID();
1015 NumProc = Comm.NumProc();
1016
1017 NumMyRows = A.NumMyRows();
1018 //NumMyCols = A.NumMyCols();
1019
1020 NumGlobalRows = A.NumGlobalRows64();
1021 NumGlobalCols = A.NumGlobalCols64();
1022
1023 if (NumGlobalRows != NumGlobalCols)
1024 IFPACK_CHK_ERR(-1); // never tested
1025
1026 /* to be changed for rect matrices */
1027 maxdim = (NumGlobalRows>NumGlobalCols)?NumGlobalRows:NumGlobalCols;
1028 maxdim /= NumPDEEqns;
1029
1030 m = 1 + maxdim;
1031 nr = NumGlobalRows / NumPDEEqns + 1;
1032 nc = NumGlobalCols / NumPDEEqns + 1;
1033
1034 if (munt == 'E') {
1035 u2dot = 72.0/conv;
1036 paperx = 21.0;
1037 siz = 10.0;
1038 }
1039 else {
1040 u2dot = 72.0;
1041 paperx = 8.5*conv;
1042 siz = siz*conv;
1043 }
1044
1045 /* left and right margins (drawing is centered) */
1046
1047 lrmrgn = (paperx-siz)/2.0;
1048
1049 /* bottom margin : 2 cm */
1050
1051 botmrgn = 2.0;
1052 /* c scaling factor */
1053 scfct = siz*u2dot/m;
1054 /* matrix frame line witdh */
1055 frlw = 0.25;
1056 /* font size for title (cm) */
1057 fnstit = 0.5;
1058 /* mfh 23 Jan 2013: title is always nonnull, since it's an array of
1059 fixed nonzero length. The 'if' test thus results in a compiler
1060 warning. */
1061 /*if (title) ltit = strlen(title);*/
1062 /*else ltit = 0;*/
1063 ltit = strlen(title);
1064
1065 /* position of title : centered horizontally */
1066 /* at 1.0 cm vertically over the drawing */
1067 ytitof = 1.0;
1068 xtit = paperx/2.0;
1069 ytit = botmrgn+siz*nr/m + ytitof;
1070 /* almost exact bounding box */
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;
1075 if (ltit == 0) {
1076 yt = yt + (ytitof+fnstit*0.70)*u2dot;
1077 }
1078 /* add some room to bounding box */
1079 delt = 10.0;
1080 xl = xl-delt;
1081 xr = xr+delt;
1082 yb = yb-delt;
1083 yt = yt+delt;
1084
1085 /* correction for title under the drawing */
1086 if ((ptitle == 0) && (ltit == 0)) {
1087 ytit = botmrgn + fnstit*0.3;
1088 botmrgn = botmrgn + ytitof + fnstit*0.7;
1089 }
1090
1091 /* begin of output */
1092
1093 if (MyPID == 0) {
1094
1095 fp = fopen(FileName,"w");
1096
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",
1100 xl,yb,xr,yt);
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");
1107
1108 /* we leave margins etc. in cm so it is easy to modify them if
1109 needed by editing the output file */
1110 fprintf(fp,"%s","gsave\n");
1111 if (ltit != 0) {
1112 fprintf(fp,"/Helvetica findfont %e cm scalefont setfont\n",
1113 fnstit);
1114 fprintf(fp,"%f cm %f cm moveto\n",
1115 xtit,ytit);
1116 fprintf(fp,"(%s) Cshow\n", title);
1117 fprintf(fp,"%f cm %f cm translate\n",
1118 lrmrgn,botmrgn);
1119 }
1120 fprintf(fp,"%f cm %d div dup scale \n",
1121 siz, (int) m);
1122 /* draw a frame around the matrix */
1123
1124 fprintf(fp,"%f setlinewidth\n",
1125 frlw);
1126 fprintf(fp,"%s","newpath\n");
1127 fprintf(fp,"%s","0 0 moveto ");
1128 if (square) {
1129 printf("------------------- %d\n", (int) m);
1130 fprintf(fp,"%d %d lineto\n",
1131 (int) m, 0);
1132 fprintf(fp,"%d %d lineto\n",
1133 (int) m, (int) m);
1134 fprintf(fp,"%d %d lineto\n",
1135 0, (int) m);
1136 }
1137 else {
1138 fprintf(fp,"%d %d lineto\n",
1139 (int) nc, 0);
1140 fprintf(fp,"%d %d lineto\n",
1141 (int) nc, (int) nr);
1142 fprintf(fp,"%d %d lineto\n",
1143 0, (int) nr);
1144 }
1145 fprintf(fp,"%s","closepath stroke\n");
1146
1147 /* plotting loop */
1148
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");
1153
1154 fclose(fp);
1155 }
1156
1157 int MaxEntries = A.MaxNumEntries();
1158 std::vector<int> Indices(MaxEntries);
1159 std::vector<double> Values(MaxEntries);
1160
1161 for (int pid = 0 ; pid < NumProc ; ++pid) {
1162
1163 if (pid == MyPID) {
1164
1165 fp = fopen(FileName,"a");
1166 TEUCHOS_ASSERT(fp != NULL);
1167
1168 for (int i = 0 ; i < NumMyRows ; ++i) {
1169
1170 if (i % NumPDEEqns) continue;
1171
1172 int Nnz;
1173 A.ExtractMyRowCopy(i,MaxEntries,Nnz,&Values[0],&Indices[0]);
1174
1175 long long grow = A.RowMatrixRowMap().GID64(i);
1176
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]);
1181 grow /= NumPDEEqns;
1182 gcol /= NumPDEEqns;
1183 fprintf(fp,"%lld %lld p\n",
1184 gcol, NumGlobalRows - grow - 1);
1185 }
1186 }
1187 }
1188
1189 fprintf(fp,"%s","%end of data for this process\n");
1190
1191 if( pid == NumProc - 1 )
1192 fprintf(fp,"%s","showpage\n");
1193
1194 fclose(fp);
1195 }
1196 Comm.Barrier();
1197 }
1198
1199 return(0);
1200}
Insert
Copy
int MyGlobalElements(int *MyGlobalElementList) const
int GID(int LID) const
bool GlobalIndicesInt() const
int NumMyElements() const
bool GlobalIndicesLongLong() const
virtual int MaxAll(double *PartialMaxs, double *GlobalMaxs, int Count) const=0
virtual int NumProc() const=0
virtual int MinAll(double *PartialMins, double *GlobalMins, int Count) const=0
virtual int MyPID() const=0
virtual void Barrier() const=0
virtual int SumAll(double *PartialSums, double *GlobalSums, int Count) const=0
const Epetra_Import * Importer() const
int FillComplete(bool OptimizeDataStorage=true)
double NormOne() const
int Scale(double ScalarConstant)
double NormInf() const
int Import(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
const Epetra_Comm & Comm() const
int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
int SumIntoGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
const Epetra_BlockMap & TargetMap() const
int NumVectors() const
int Abs(const Epetra_MultiVector &A)
int MinValue(double *Result) const
int MyLength() const
int MeanValue(double *Result) const
int MaxValue(double *Result) const
long long GlobalLength64() const
int PutScalar(double ScalarConstant)
virtual const char * Label() const
virtual const Epetra_Comm & Comm() const=0
virtual const char * Label() const=0
virtual const Epetra_Map & OperatorDomainMap() const=0
virtual const Epetra_Map & OperatorRangeMap() const=0
virtual int NumMyRows() const=0
virtual int Multiply(bool TransA, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const=0
virtual const Epetra_Map & RowMatrixColMap() const=0
virtual const Epetra_Map & RowMatrixRowMap() const=0
virtual int MaxNumEntries() const=0
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const=0
virtual const Epetra_Import * RowMatrixImporter() const=0
virtual double NormOne() const=0
virtual double NormInf() const=0