5 #include "Epetra_MpiComm.h" 7 #include "Epetra_SerialComm.h" 9 #include "Epetra_Map.h" 10 #include "Epetra_Vector.h" 11 #include "Epetra_Time.h" 12 #include "Epetra_Util.h" 13 #include "Epetra_CrsMatrix.h" 16 #include "Teuchos_ParameterList.hpp" 20 const Epetra_RowMatrix& A,
21 const Epetra_MultiVector& x,
22 const Epetra_MultiVector& b,
23 const Epetra_MultiVector& x_exact)
25 std::vector<double> Norm;
26 int NumVectors = x.NumVectors();
27 Norm.resize(NumVectors);
28 Epetra_MultiVector Ax(x);
29 A.Multiply(
false,x,Ax);
30 Ax.Update(1.0,b,-1.0);
32 bool TestPassed =
false;
33 double TotalNorm = 0.0;
34 for (
int i = 0 ; i < NumVectors ; ++i) {
37 if (
verbose && A.Comm().MyPID() == 0)
38 std::cout <<
"||Ax - b|| = " << TotalNorm << std::endl;
39 if (TotalNorm < 1e-5 )
44 Ax.Update (1.0,x,-1.0,x_exact,0.0);
46 for (
int i = 0 ; i < NumVectors ; ++i) {
49 if (
verbose && A.Comm().MyPID() == 0)
50 std::cout <<
"||Ax - b|| = " << TotalNorm << std::endl;
51 if (TotalNorm < 1e-5 )
61 int main(
int argc,
char *argv[]) {
64 MPI_Init(&argc, &argv);
65 Epetra_MpiComm Comm(MPI_COMM_WORLD);
67 Epetra_SerialComm Comm;
70 int NumGlobalElements = 1000;
75 if ( argc > 1 && argv[1][0] ==
'-' && argv[1][1] ==
'q' )
verbose = false ;
81 int* part =
new int[NumGlobalElements];
83 if (Comm.MyPID() == 0) {
86 for(
int i=0 ; i<NumGlobalElements ; ++i ) {
87 unsigned int r = Util.RandomInt();
88 part[i] = r%(Comm.NumProc());
92 Comm.Broadcast(part,NumGlobalElements,0);
95 int NumMyElements = 0;
96 for (
int i = 0 ; i < NumGlobalElements ; ++i) {
97 if (part[i] == Comm.MyPID())
102 int* MyGlobalElements =
new int[NumMyElements];
104 for (
int i = 0 ; i < NumGlobalElements ; ++i) {
105 if (part[i] == Comm.MyPID() )
106 MyGlobalElements[count++] = i;
109 Epetra_Map Map(NumGlobalElements,NumMyElements,MyGlobalElements,
118 Epetra_CrsMatrix Matrix(Copy,Map,NumGlobalElements);
120 int* Indices =
new int[NumGlobalElements];
121 double* Values =
new double[NumGlobalElements];
123 for (
int i = 0 ; i < NumGlobalElements ; ++i)
126 for (
int i = 0 ; i < NumMyElements ; ++i)
128 int iGlobal = MyGlobalElements[i];
129 const int MakeNotDense = 1;
130 int Min_jGlobal = std::min(i,MakeNotDense );
131 for (
int jGlobal = Min_jGlobal ; jGlobal < NumGlobalElements ; ++jGlobal) {
132 if (iGlobal == jGlobal)
133 Values[jGlobal-Min_jGlobal] = 1.0 * (NumGlobalElements + 1 ) *
134 (NumGlobalElements + 1);
135 else if (iGlobal > jGlobal)
136 Values[jGlobal-Min_jGlobal] = -1.0*(jGlobal+1);
138 Values[jGlobal-Min_jGlobal] = 1.0*(iGlobal+1);
140 Matrix.InsertGlobalValues(MyGlobalElements[i],
141 NumGlobalElements-MakeNotDense, Values, &Indices[Min_jGlobal]);
145 Matrix.FillComplete();
147 delete[] MyGlobalElements;
156 Epetra_MultiVector x(Map,NumVectors);
157 Epetra_MultiVector x_exact(Map,NumVectors);
158 Epetra_MultiVector b(Map,NumVectors);
166 Epetra_LinearProblem Problem(&Matrix, &x, &b);
169 Teuchos::ParameterList List;
170 List.set(
"MaxProcs", Comm.NumProc());
172 Solver->SetParameters(List);
186 return(EXIT_SUCCESS);
Amesos_Superludist: An object-oriented wrapper for Superludist.
bool CheckError(bool verbose, const Epetra_RowMatrix &A, const Epetra_MultiVector &x, const Epetra_MultiVector &b, const Epetra_MultiVector &x_exact)
Amesos_TestRowMatrix: a class to test Epetra_RowMatrix based codes.
#define AMESOS_CHK_ERR(a)
virtual int Multiply(bool TransA, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of a Epetra_RowMatrix multiplied by a Epetra_MultiVector X in Y.
int main(int argc, char *argv[])