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