Amesos Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
Test_DSCPACK/cxx_main.cpp
Go to the documentation of this file.
1//
2// This test does not exercise multiple processes on DSCPACK.
3// DSCPACK uses only one process for dense matrices.
4//
5//
6#include "Amesos_ConfigDefs.h"
7
8#ifdef HAVE_MPI
9#include "mpi.h"
10#include "Epetra_MpiComm.h"
11#else
12#include "Epetra_SerialComm.h"
13#endif
14#include "Epetra_Map.h"
15#include "Epetra_Vector.h"
16#include "Epetra_Time.h"
17#include "Epetra_Util.h"
18#include "Epetra_CrsMatrix.h"
19#include "Amesos_Dscpack.h"
21#include "Teuchos_ParameterList.hpp"
22
23//=============================================================================
24bool CheckError(const Epetra_RowMatrix& A,
25 const Epetra_MultiVector& x,
26 const Epetra_MultiVector& b,
27 const Epetra_MultiVector& x_exact)
28{
29 std::vector<double> Norm;
30 int NumVectors = x.NumVectors();
31 Norm.resize(NumVectors);
32 Epetra_MultiVector Ax(x);
33 A.Multiply(false,x,Ax);
34 Ax.Update(1.0,b,-1.0);
35 Ax.Norm2(&Norm[0]);
36 bool TestPassed = false;
37 double TotalNorm = 0.0;
38 for (int i = 0 ; i < NumVectors ; ++i) {
39 TotalNorm += Norm[i];
40 }
41 if (A.Comm().MyPID() == 0)
42 std::cout << "||Ax - b|| = " << TotalNorm << std::endl;
43 if (TotalNorm < 1e-5 )
44 TestPassed = true;
45 else
46 TestPassed = false;
47
48 Ax.Update (1.0,x,-1.0,x_exact,0.0);
49 Ax.Norm2(&Norm[0]);
50 for (int i = 0 ; i < NumVectors ; ++i) {
51 TotalNorm += Norm[i];
52 }
53 if (A.Comm().MyPID() == 0)
54 std::cout << "||Ax - b|| = " << TotalNorm << std::endl;
55 if (TotalNorm < 1e-5 )
56 TestPassed = true;
57 else
58 TestPassed = false;
59
60 return(TestPassed);
61}
62
63//=============================================================================
64int main(int argc, char *argv[]) {
65
66#ifdef HAVE_MPI
67 MPI_Init(&argc, &argv);
68 Epetra_MpiComm Comm(MPI_COMM_WORLD);
69#else
70 Epetra_SerialComm Comm;
71#endif
72
73 int NumGlobalElements = 1000;
74 int NumVectors = 7;
75
76 // =================== //
77 // create a random map //
78 // =================== //
79
80 int* part = new int[NumGlobalElements];
81
82 if (Comm.MyPID() == 0) {
83 Epetra_Util Util;
84
85 for( int i=0 ; i<NumGlobalElements ; ++i ) {
86 unsigned int r = Util.RandomInt();
87 part[i] = r%(Comm.NumProc());
88 }
89 }
90
91 Comm.Broadcast(part,NumGlobalElements,0);
92
93 // count the elements assigned to this proc
94 int NumMyElements = 0;
95 for (int i = 0 ; i < NumGlobalElements ; ++i) {
96 if (part[i] == Comm.MyPID())
97 NumMyElements++;
98 }
99
100 // get the loc2global list
101 int* MyGlobalElements = new int[NumMyElements];
102 int count = 0;
103 for (int i = 0 ; i < NumGlobalElements ; ++i) {
104 if (part[i] == Comm.MyPID() )
105 MyGlobalElements[count++] = i;
106 }
107
108 Epetra_Map Map(NumGlobalElements,NumMyElements,MyGlobalElements,
109 0,Comm);
110
111 delete [] part;
112
113 // ===================== //
114 // Create a dense matrix //
115 // ===================== //
116
117 Epetra_CrsMatrix Matrix(Copy,Map,NumGlobalElements);
118
119 int* Indices = new int[NumGlobalElements];
120 double* Values = new double[NumGlobalElements];
121
122 for (int i = 0 ; i < NumGlobalElements ; ++i)
123 Indices[i] = i;
124
125 for (int i = 0 ; i < NumMyElements ; ++i) {
126 int iGlobal = MyGlobalElements[i];
127 for (int jGlobal = 0 ; jGlobal < NumGlobalElements ; ++jGlobal) {
128 if (iGlobal == jGlobal)
129 Values[jGlobal] = 1.0 * (NumGlobalElements + 1 ) *
130 (NumGlobalElements + 1);
131 else if (iGlobal > jGlobal)
132 Values[jGlobal] = 1.0*(jGlobal+1);
133 else
134 Values[jGlobal] = 1.0*(iGlobal+1);
135 }
136 Matrix.InsertGlobalValues(MyGlobalElements[i],
137 NumGlobalElements, Values, Indices);
138
139 }
140
141 Matrix.FillComplete();
142
143 delete [] MyGlobalElements;
144 delete [] Indices;
145 delete [] Values;
146
147 // ======================== //
148 // other data for this test //
149 // ======================== //
150
151 Amesos_TestRowMatrix A(&Matrix);
152 Epetra_MultiVector x(Map,NumVectors);
153 Epetra_MultiVector x_exact(Map,NumVectors);
154 Epetra_MultiVector b(Map,NumVectors);
155 x_exact.Random();
156 Matrix.Multiply(false,x_exact,b);
157
158 // =========== //
159 // AMESOS PART //
160 // =========== //
161
162 Epetra_LinearProblem Problem;
163 Amesos_Dscpack Solver(Problem);
164
165 Problem.SetOperator(&A);
166 Problem.SetLHS(&x);
167 Problem.SetRHS(&b);
168
169 Teuchos::ParameterList ParamList ;
170 ParamList.set("MaxProcs", Comm.NumProc());
171 EPETRA_CHK_ERR(Solver.SetParameters(ParamList));
172
175 AMESOS_CHK_ERR(Solver.Solve());
176
177 if (!CheckError(Matrix,x,b,x_exact))
178 AMESOS_CHK_ERR(-1);
179
180#ifdef HAVE_MPI
181 MPI_Finalize();
182#endif
183
184 return 0;
185}
#define AMESOS_CHK_ERR(a)
#define EPETRA_CHK_ERR(xxx)
bool CheckError(const Epetra_RowMatrix &A, const Epetra_MultiVector &x, const Epetra_MultiVector &b, const Epetra_MultiVector &x_exact)
int main(int argc, char *argv[])
Amesos_Dscpack: An object-oriented wrapper for Dscpack.
int NumericFactorization()
Performs NumericFactorization on the matrix A.
int Solve()
Solves A X = B (or AT x = B)
int SetParameters(Teuchos::ParameterList &ParameterList)
Updates internal variables.
int SymbolicFactorization()
Performs SymbolicFactorization on the matrix A.
Amesos_TestRowMatrix: a class to test Epetra_RowMatrix based codes.