Ifpack Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
test/TestAll/cxx_main.cpp
Go to the documentation of this file.
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
45#ifdef HAVE_MPI
46#include "Epetra_MpiComm.h"
47#else
48#include "Epetra_SerialComm.h"
49#endif
50#include "Epetra_CrsMatrix.h"
51#include "Epetra_Vector.h"
53#include "Galeri_Maps.h"
54#include "Galeri_CrsMatrices.h"
55#include "Teuchos_ParameterList.hpp"
56#include "Teuchos_RefCountPtr.hpp"
58#include "AztecOO.h"
63#include "Ifpack_Amesos.h"
64#include "Ifpack_Utils.h"
65#include "Ifpack_Chebyshev.h"
66#include "Ifpack_Polynomial.h"
67#include "Ifpack_Krylov.h"
68
69template <class T>
70bool Test(const Teuchos::RefCountPtr<Epetra_RowMatrix>& Matrix, Teuchos::ParameterList& List)
71{
72
73 int NumVectors = 1;
74 bool UseTranspose = false;
75
76 Epetra_MultiVector LHS(Matrix->OperatorDomainMap(),NumVectors);
77 Epetra_MultiVector RHS(Matrix->OperatorRangeMap(),NumVectors);
78 Epetra_MultiVector LHSexact(Matrix->OperatorDomainMap(),NumVectors);
79
80 LHS.PutScalar(0.0);
81 LHSexact.Random();
82 Matrix->Multiply(UseTranspose,LHSexact,RHS);
83
84 Epetra_LinearProblem Problem(&*Matrix,&LHS,&RHS);
85
86 Teuchos::RefCountPtr<T> Prec;
87
88 Prec = Teuchos::rcp( new T(&*Matrix) );
89 assert(Prec != Teuchos::null);
90
91 IFPACK_CHK_ERR(Prec->SetParameters(List));
92 IFPACK_CHK_ERR(Prec->Initialize());
93 IFPACK_CHK_ERR(Prec->Compute());
94
95 // create the AztecOO solver
96 AztecOO AztecOOSolver(Problem);
97
98 // specify solver
99 AztecOOSolver.SetAztecOption(AZ_solver,AZ_gmres);
100 AztecOOSolver.SetAztecOption(AZ_output,32);
101
102 AztecOOSolver.SetPrecOperator(&*Prec);
103
104 // solver. The solver should converge in one iteration,
105 // or maximum two (numerical errors)
106 AztecOOSolver.Iterate(1550,1e-8);
107
108 cout << *Prec;
109
110 std::vector<double> Norm(NumVectors);
111 LHS.Update(1.0,LHSexact,-1.0);
112 LHS.Norm2(&Norm[0]);
113 for (int i = 0 ; i < NumVectors ; ++i) {
114 cout << "Norm[" << i << "] = " << Norm[i] << endl;
115 if (Norm[i] > 1e-3)
116 return(false);
117 }
118 return(true);
119
120}
121
122int main(int argc, char *argv[])
123{
124
125#ifdef HAVE_MPI
126 MPI_Init(&argc,&argv);
127 Epetra_MpiComm Comm(MPI_COMM_WORLD);
128#else
130#endif
131
132 bool verbose = (Comm.MyPID() == 0);
133
134 Teuchos::ParameterList GaleriList;
135 const int nx = 30;
136 GaleriList.set("n", nx * nx);
137 GaleriList.set("nx", nx);
138 GaleriList.set("ny", nx);
139 Teuchos::RefCountPtr<Epetra_Map> Map = Teuchos::rcp( Galeri::CreateMap("Linear", Comm, GaleriList) );
140 Teuchos::RefCountPtr<Epetra_RowMatrix> Matrix = Teuchos::rcp( Galeri::CreateCrsMatrix("Laplace2D", &*Map, GaleriList) );
141
142 Teuchos::ParameterList List, DefaultList;
143
144 // test the preconditioner
145 int TestPassed = true;
146
147 if (!Test<Ifpack_Chebyshev>(Matrix,List))
148 {
149 TestPassed = false;
150 }
151
152 List.set("polynomial: degree",3);
153 if (!Test<Ifpack_Polynomial>(Matrix,List))
154 {
155 TestPassed = false;
156 }
157
158 List = DefaultList;
159 List.set("krylov: tolerance", 1e-14);
160 List.set("krylov: iterations", 100);
161 List.set("krylov: preconditioner", 2);
162 List.set("krylov: block size", 9);
163 List.set("krylov: number of sweeps", 2);
164 if (!Test<Ifpack_Krylov>(Matrix,List))
165 {
166 TestPassed = false;
167 }
168
169 if (!Test< Ifpack_AdditiveSchwarz<Ifpack_Krylov> >(Matrix,List)) {
170 TestPassed = false;
171 }
172
173
174 if (!Test<Ifpack_Amesos>(Matrix,List))
175 {
176 TestPassed = false;
177 }
178
179 // FIXME
180#if 0
181 if (!Test<Ifpack_AdditiveSchwarz<Ifpack_BlockRelaxation<Ifpack_DenseContainer> > >(Matrix,List)) {
182 TestPassed = false;
183 }
184#endif
185
186 // this is ok as long as just one sweep is applied
187 List = DefaultList;
188 List.set("relaxation: type", "Gauss-Seidel");
189 if (!Test<Ifpack_PointRelaxation>(Matrix,List)) {
190 TestPassed = false;
191 }
192
193 // this is ok as long as just one sweep is applied
194 List = DefaultList;
195 List.set("relaxation: type", "symmetric Gauss-Seidel");
196 List.set("relaxation: sweeps", 5);
197 List.set("partitioner: local parts", 128);
198 List.set("partitioner: type", "linear");
199 if (!Test<Ifpack_BlockRelaxation<Ifpack_DenseContainer> >(Matrix,List)) {
200 TestPassed = false;
201 }
202
203 // this is ok as long as just one sweep is applied
204 List = DefaultList;
205 List.set("relaxation: type", "symmetric Gauss-Seidel");
206 List.set("partitioner: local parts", 128);
207 List.set("partitioner: type", "linear");
208 if (!Test<Ifpack_BlockRelaxation<Ifpack_SparseContainer<Ifpack_Amesos> > >(Matrix,List)) {
209 TestPassed = false;
210 }
211
212 // this is ok as long as just one sweep is applied
213 List = DefaultList;
214 List.set("relaxation: type", "symmetric Gauss-Seidel");
215 List.set("partitioner: local parts", 128);
216 List.set("partitioner: type", "linear");
217 if (!Test<Ifpack_AdditiveSchwarz<Ifpack_BlockRelaxation<Ifpack_SparseContainer<Ifpack_Amesos> > > >(Matrix,List)) {
218 TestPassed = false;
219 }
220 if (!TestPassed) {
221 cerr << "Test `TestAll.exe' FAILED!" << endl;
222 exit(EXIT_FAILURE);
223 }
224
225#ifdef HAVE_MPI
226 MPI_Finalize();
227#endif
228 if (verbose)
229 cout << "Test `TestAll.exe' passed!" << endl;
230
231 return(EXIT_SUCCESS);
232}
#define IFPACK_CHK_ERR(ifpack_err)
#define RHS(a)
Definition MatGenFD.c:60
int MyPID() const
int Update(double ScalarA, const Epetra_MultiVector &A, double ScalarThis)
int Norm2(double *Result) const
int PutScalar(double ScalarConstant)
Ifpack_AdditiveSchwarz: a class to define Additive Schwarz preconditioners of Epetra_RowMatrix's.
Ifpack_BlockRelaxation: a class to define block relaxation preconditioners of Epetra_RowMatrix's.
const int NumVectors
int main(int argc, char *argv[])
bool Test(const Teuchos::RefCountPtr< Epetra_RowMatrix > &Matrix, Teuchos::ParameterList &List)