Amesos Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
Test_Singular/cxx_main.cpp
Go to the documentation of this file.
1#include "Amesos_ConfigDefs.h"
2
3#ifdef HAVE_MPI
4#include "mpi.h"
5#include "Epetra_MpiComm.h"
6#else
7#include "Epetra_SerialComm.h"
8#endif
9#include "Epetra_Util.h"
10#include "Epetra_Map.h"
11#include "Epetra_Vector.h"
12#include "Epetra_CrsMatrix.h"
13#include "Epetra_LinearProblem.h"
14#include "Epetra_Export.h"
15#include "Amesos.h"
16#include "Teuchos_ParameterList.hpp"
17#include "Teuchos_RCP.hpp"
18#include "Galeri_ReadHB.h"
19
20//============ //
21// main driver //
22//============ //
23
24int main(int argc, char *argv[])
25{
26#ifdef HAVE_MPI
27 MPI_Init(&argc, &argv);
28 Epetra_MpiComm Comm(MPI_COMM_WORLD);
29#else
30 Epetra_SerialComm Comm;
31#endif
32
33// ================= //
34 // reading HB matrix //
35 // ================= //
36
37 // HB files are for serial matrices. Hence, only
38 // process 0 reads this matrix (and if present
39 // solution and RHS). Then, this matrix will be redistributed
40 // using epetra capabilities.
41 // All variables that begin with "read" refer to the
42 // HB matrix read by process 0.
43 Epetra_Map* readMap;
44 Epetra_CrsMatrix* readA;
45 Epetra_Vector* readx;
46 Epetra_Vector* readb;
47 Epetra_Vector* readxexact;
48
49 // Name of input file with a numerical singularity
50 std::string matrix_file="bcsstm05_ns.rua";
51 try
52 {
53 Galeri::ReadHB(matrix_file.c_str(), Comm, readMap,
54 readA, readx, readb, readxexact);
55 }
56 catch(...)
57 {
58 std::cout << "Caught exception, maybe file name is incorrect" << std::endl;
59#ifdef HAVE_MPI
60 MPI_Finalize();
61#else
62 // not to break test harness
63 exit(EXIT_SUCCESS);
64#endif
65 }
66
67 // Create uniform distributed map.
68 // Note that linear map are used for simplicity only!
69 // Amesos (through Epetra) can support *any* map.
70
71 Epetra_Map* mapPtr = 0;
72 if(readMap->GlobalIndicesInt())
73 mapPtr = new Epetra_Map((int) readMap->NumGlobalElements(), 0, Comm);
74 else if(readMap->GlobalIndicesLongLong())
75 mapPtr = new Epetra_Map(readMap->NumGlobalElements(), 0, Comm);
76 else
77 assert(false);
78
79 Epetra_Map& map = *mapPtr;
80
81 // Create the distributed matrix, based on Map.
82 Epetra_CrsMatrix A(Copy, map, 0);
83
84 const Epetra_Map &OriginalMap = readA->RowMatrixRowMap() ;
85 assert (OriginalMap.SameAs(*readMap));
86 Epetra_Export exporter(OriginalMap, map);
87
88 Epetra_Vector x(map); // distributed solution
89 Epetra_Vector b(map); // distributed rhs
90 Epetra_Vector xexact(map); // distributed exact solution
91
92 // Exports from data defined on processor 0 to distributed.
93 x.Export(*readx, exporter, Add);
94 b.Export(*readb, exporter, Add);
95 xexact.Export(*readxexact, exporter, Add);
96 A.Export(*readA, exporter, Add);
97 A.FillComplete();
98
99 // deletes memory
100 delete readMap;
101 delete readA;
102 delete readx;
103 delete readb;
104 delete readxexact;
105
106 // Creates an epetra linear problem, contaning matrix
107 // A, solution x and rhs b.
108 Epetra_LinearProblem problem(&A,&x,&b);
109
110 // =========== //
111 // AMESOS PART //
112 // =========== //
113
114 // Create the factory
115 Amesos factory;
116
117 // Create the solver
118 std::string solverType = "Klu";
119 Teuchos::RCP<Amesos_BaseSolver> solver = Teuchos::rcp( factory.Create(solverType, problem) );
120
121 // Perform symbolic factorization
122 int symRet = solver->SymbolicFactorization();
123 if (symRet != 0) {
124 std::cout << "Processor "<< map.Comm().MyPID() << " : Symbolic factorization did not complete!" << std::endl;
125 }
126
127 // Perform numeric factorization
128 int numRet = solver->NumericFactorization();
129 if (numRet != 0) {
130 std::cout << "Processor "<< map.Comm().MyPID() << " : Numeric factorization did not complete!" << std::endl;
131 }
132
133 // Check that all processors returned error -22 (Numerically Singular)!
134 int minRet = 0;
135 Comm.MinAll( &numRet, &minRet, 1 );
136
137 if (minRet != NumericallySingularMatrixError) {
138 if (Comm.MyPID()==0)
139 std::cout << std::endl << "End Result: TEST FAILED" << std::endl;
140 }
141 else {
142 if (Comm.MyPID()==0)
143 std::cout << std::endl << "End Result: TEST PASSED" << std::endl;
144 }
145
146 delete mapPtr;
147
148#ifdef HAVE_MPI
149 MPI_Finalize();
150#endif
151
152 return(0);
153}
const int NumericallySingularMatrixError
int main(int argc, char *argv[])
Factory for binding a third party direct solver to an Epetra_LinearProblem.
Definition Amesos.h:44