Ifpack Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
test/IHSS-SORa/cxx_main.cpp
Go to the documentation of this file.
1
2/*@HEADER
3// ***********************************************************************
4//
5// Ifpack: Object-Oriented Algebraic Preconditioner Package
6// Copyright (2002) Sandia Corporation
7//
8// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
9// license for use of this work by or on behalf of the U.S. Government.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
39//
40// ***********************************************************************
41//@HEADER
42*/
43
44#include "Ifpack_ConfigDefs.h"
45
46#include "Epetra_ConfigDefs.h"
47#ifdef HAVE_MPI
48#include "mpi.h"
49#include "Epetra_MpiComm.h"
50#else
51#include "Epetra_SerialComm.h"
52#endif
53#include "Epetra_Comm.h"
54#include "Epetra_Map.h"
55#include "Epetra_Time.h"
56#include "Epetra_BlockMap.h"
57#include "Epetra_MultiVector.h"
58#include "Epetra_Vector.h"
59#include "Epetra_Export.h"
60#include "AztecOO.h"
61#include "Galeri_Maps.h"
62#include "Galeri_CrsMatrices.h"
63#include "Ifpack_IHSS.h"
64#include "Ifpack_SORa.h"
65#include "Ifpack.h"
66#include "Teuchos_RefCountPtr.hpp"
67
68// function for fancy output
69
70std::string toString(const int& x) {
71 char s[100];
72 sprintf(s, "%d", x);
73 return std::string(s);
74}
75
76std::string toString(const double& x) {
77 char s[100];
78 sprintf(s, "%g", x);
79 return std::string(s);
80}
81
82// main driver
83
84int main(int argc, char *argv[]) {
85
86#ifdef HAVE_MPI
87 MPI_Init(&argc,&argv);
88 Epetra_MpiComm Comm (MPI_COMM_WORLD);
89#else
91#endif
92
93 Teuchos::ParameterList GaleriList;
94 int nx = 30;
95
96 GaleriList.set("nx", nx);
97 // GaleriList.set("ny", nx * Comm.NumProc());
98 GaleriList.set("ny", nx);
99 GaleriList.set("mx", 1);
100 GaleriList.set("my", Comm.NumProc());
101 GaleriList.set("alpha", .0);
102 GaleriList.set("diff", 1.0);
103 GaleriList.set("conv", 100.0);
104
105 Teuchos::RefCountPtr<Epetra_Map> Map = Teuchos::rcp( Galeri::CreateMap("Cartesian2D", Comm, GaleriList) );
106 Teuchos::RefCountPtr<Epetra_CrsMatrix> A = Teuchos::rcp( Galeri::CreateCrsMatrix("UniFlow2D", &*Map, GaleriList) );
107 Teuchos::RefCountPtr<Epetra_MultiVector> LHS = Teuchos::rcp( new Epetra_MultiVector(*Map, 1) );
108 Teuchos::RefCountPtr<Epetra_MultiVector> RHS = Teuchos::rcp( new Epetra_MultiVector(*Map, 1) );
109 LHS->PutScalar(0.0); RHS->Random();
110 Ifpack Factory;
111 int Niters = 100;
112
113 // ============================= //
114 // Construct IHSS preconditioner //
115 // ============================= //
116 Teuchos::RefCountPtr<Ifpack_Preconditioner> Prec = Teuchos::rcp( Factory.Create("IHSS", &*A,0) );
117 Teuchos::ParameterList List;
118 List.set("ihss: hermetian type","ILU");
119 List.set("ihss: skew hermetian type","ILU");
120 List.set("ihss: ratio eigenvalue",100.0);
121 // Could set sublist values here to better control the ILU, but this isn't needed for this example.
122 IFPACK_CHK_ERR(Prec->SetParameters(List));
123 IFPACK_CHK_ERR(Prec->Compute());
124
125 // ============================= //
126 // Create solver Object //
127 // ============================= //
128
129 AztecOO solver;
130 solver.SetUserMatrix(&*A);
131 solver.SetLHS(&*LHS);
132 solver.SetRHS(&*RHS);
133 solver.SetAztecOption(AZ_solver,AZ_gmres);
134 solver.SetPrecOperator(&*Prec);
135 solver.SetAztecOption(AZ_output, 1);
136 solver.Iterate(Niters, 1e-8);
137
138 // ============================= //
139 // Construct SORa preconditioner //
140 // ============================= //
141 Teuchos::RefCountPtr<Ifpack_Preconditioner> Prec2 = Teuchos::rcp( Factory.Create("SORa", &*A,0) );
142 Teuchos::ParameterList List2;
143 List2.set("sora: sweeps",1);
144 List2.set("sora: use global damping",true);
145 List2.set("sora: eigen-analysis: random seed",(unsigned int)24601);
146 // Could set sublist values here to better control the ILU, but this isn't needed for this example.
147 IFPACK_CHK_ERR(Prec2->SetParameters(List2));
148 IFPACK_CHK_ERR(Prec2->Compute());
149
150 // ============================= //
151 // Create solver Object //
152 // ============================= //
153 AztecOO solver2;
154 LHS->PutScalar(0.0);
155 solver2.SetUserMatrix(&*A);
156 solver2.SetLHS(&*LHS);
157 solver2.SetRHS(&*RHS);
158 solver2.SetAztecOption(AZ_solver,AZ_gmres);
159 solver2.SetPrecOperator(&*Prec2);
160 solver2.SetAztecOption(AZ_output, 1);
161 solver2.Iterate(Niters, 1e-8);
162
163 // ============================= //
164 // Construct a second SORa preconditioner to check seeds //
165 // ============================= //
166 Teuchos::RefCountPtr<Ifpack_Preconditioner> Prec3 = Teuchos::rcp( Factory.Create("SORa", &*A,0) );
167 Teuchos::ParameterList List3;
168 List3.set("sora: sweeps",1);
169 List3.set("sora: use global damping",true);
170 List3.set("sora: eigen-analysis: random seed",(unsigned int)24601);
171 // Could set sublist values here to better control the ILU, but this isn't needed for this example.
172 IFPACK_CHK_ERR(Prec3->SetParameters(List2));
173 IFPACK_CHK_ERR(Prec3->Compute());
174
175 Teuchos::RCP<Ifpack_SORa> Prec2_SORa = Teuchos::rcp_dynamic_cast<Ifpack_SORa>(Prec2);
176 Teuchos::RCP<Ifpack_SORa> Prec3_SORa = Teuchos::rcp_dynamic_cast<Ifpack_SORa>(Prec3);
177 double diff = Prec2_SORa->GetLambdaMax()-Prec3_SORa->GetLambdaMax();
178 if(diff > 1e-12) return EXIT_FAILURE;
179
180
181#ifdef HAVE_MPI
182 MPI_Finalize() ;
183#endif
184
185 return(EXIT_SUCCESS);
186}
#define IFPACK_CHK_ERR(ifpack_err)
#define RHS(a)
Definition: MatGenFD.c:60
int NumProc() const
static Ifpack_Preconditioner * Create(EPrecType PrecType, Epetra_RowMatrix *Matrix, const int overlap=0, bool overrideSerialDefault=false)
Creates an instance of Ifpack_Preconditioner given the enum value of the preconditioner type (can not...
Definition: Ifpack.cpp:253
int main(int argc, char *argv[])
std::string toString(const int &x)