IFPACK Development
Loading...
Searching...
No Matches
Ifpack_Euclid.cpp
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_Euclid.h"
44#if defined(HAVE_EUCLID) && defined(HAVE_MPI)
45
46#include "Ifpack_Utils.h"
47#include <algorithm>
48#include "Epetra_MpiComm.h"
49#include "Epetra_IntVector.h"
50#include "Epetra_Import.h"
51#include "Teuchos_ParameterList.hpp"
52#include "Teuchos_RCP.hpp"
53#include "getRow_dh.h"
54
55using Teuchos::RCP;
56using Teuchos::rcp;
57
58Ifpack_Euclid::Ifpack_Euclid(Epetra_CrsMatrix* A):
59 A_(rcp(A,false)),
60 UseTranspose_(false),
61 Condest_(-1),
62 IsInitialized_(false),
63 IsComputed_(false),
64 Label_(),
65 NumInitialize_(0),
66 NumCompute_(0),
67 NumApplyInverse_(0),
68 InitializeTime_(0.0),
69 ComputeTime_(0.0),
70 ApplyInverseTime_(0.0),
71 ComputeFlops_(0.0),
72 ApplyInverseFlops_(0.0),
73 Time_(A_->Comm()),
74 SetLevel_(1),
75 SetBJ_(0),
76 SetStats_(0),
77 SetMem_(0),
78 SetSparse_(0.0),
79 SetRowScale_(0),
80 SetILUT_(0.0)
81{
82 // Here we need to change the view of each row to have global indices. This is
83 // because Euclid directly extracts a row view and expects global indices.
84 for(int i = 0; i < A_->NumMyRows(); i++){
85 int *indices;
86 int len;
87 A_->Graph().ExtractMyRowView(i, len, indices);
88 for(int j = 0; j < len; j++){
89 indices[j] = A_->GCID(indices[j]);
90 }
91 }
92} //Constructor
93
94//==============================================================================
95void Ifpack_Euclid::Destroy(){
96 // Destroy the euclid solver, only if it was setup
97 if(IsComputed()){
98 Euclid_dhDestroy(eu);
99 }
100 // Delete these euclid varaiables if they were created
101 if(IsInitialized()){
102 Parser_dhDestroy(parser_dh);
103 parser_dh = NULL;
104 TimeLog_dhDestroy(tlog_dh);
105 tlog_dh = NULL;
106 Mem_dhDestroy(mem_dh);
107 mem_dh = NULL;
108 }
109 // Now that Euclid is done with the matrix, we change it back to having local indices.
110 for(int i = 0; i < A_->NumMyRows(); i++){
111 int *indices;
112 int len;
113 A_->Graph().ExtractMyRowView(i, len, indices);
114 for(int j = 0; j < len; j++){
115 indices[j] = A_->LCID(indices[j]);
116 }
117 }
118} //Destroy()
119
120//==============================================================================
121int Ifpack_Euclid::Initialize(){
122 //These are global variables in Euclid
123 comm_dh = GetMpiComm();
124 MPI_Comm_size(comm_dh, &np_dh);
125 MPI_Comm_rank(comm_dh, &myid_dh);
126 Time_.ResetStartTime();
127 if(mem_dh == NULL){
128 Mem_dhCreate(&mem_dh);
129 }
130 if (tlog_dh == NULL) {
131 TimeLog_dhCreate(&tlog_dh);
132 }
133
134 if (parser_dh == NULL) {
135 Parser_dhCreate(&parser_dh);
136 }
137 Parser_dhInit(parser_dh, 0, NULL);
138 // Create the solver, this doesn't malloc anything yet, so it's only destroyed if Compute() is called.
139 Euclid_dhCreate(&eu);
140 IsInitialized_=true;
141 NumInitialize_ = NumInitialize_ + 1;
142 InitializeTime_ = InitializeTime_ + Time_.ElapsedTime();
143 return 0;
144} //Initialize()
145
146//==============================================================================
147int Ifpack_Euclid::SetParameters(Teuchos::ParameterList& list){
148 List_ = list;
149 SetLevel_ = list.get("SetLevel", (int)1);
150 SetBJ_ = list.get("SetBJ", (int)0);
151 SetStats_ = list.get("SetStats", (int)0);
152 SetMem_ = list.get("SetMem", (int)0);
153 SetSparse_ = list.get("SetSparse", (double)0.0);
154 SetRowScale_ = list.get("SetRowScale", (int)0);
155 SetILUT_ = list.get("SetILUT", (double)0.0);
156 return 0;
157} //SetParamters()
158
159//==============================================================================
160int Ifpack_Euclid::SetParameter(std::string name, int value){
161 //Convert to lowercase (so it's case insensitive)
162 locale loc;
163 for(size_t i = 0; i < name.length(); i++){
164 name[i] = (char) tolower(name[i], loc);
165 }
166 if(name.compare("setlevel") == 0){
167 SetLevel_ = value;
168 } else if(name.compare("setbj") == 0){
169 SetBJ_ = value;
170 } else if(name.compare("setstats") == 0){
171 SetStats_ = value;
172 } else if(name.compare("setmem") == 0){
173 SetMem_ = value;
174 } else if(name.compare("setrowscale") == 0){
175 SetRowScale_ = value;
176 } else {
177 using std::cout;
178 using std::endl;
179
180 cout << "\nThe string " << name << " is not an available option." << endl;
181 IFPACK_CHK_ERR(-1);
182 }
183 return 0;
184} //SetParameter() (int)
185
186//==============================================================================
187int Ifpack_Euclid::SetParameter(std::string name, double value){
188 //Convert to lowercase (so it's case insensitive)
189 locale loc;
190 for(size_t i; i < name.length(); i++){
191 name[i] = (char) tolower(name[i], loc);
192 }
193 if(name.compare("setsparse") == 0){
194 SetSparse_ = value;
195 } else if(name.compare("setilut") == 0){
196 SetILUT_ = value;
197 } else {
198 using std::cout;
199 using std::endl;
200
201 cout << "\nThe string " << name << " is not an available option." << endl;
202 IFPACK_CHK_ERR(-1);
203 }
204 return 0;
205} //SetParameter() (double)
206
207//==============================================================================
208int Ifpack_Euclid::Compute(){
209 if(IsInitialized() == false){
210 IFPACK_CHK_ERR(Initialize());
211 }
212 Time_.ResetStartTime();
213 sprintf(Label_, "IFPACK_Euclid (level=%d, bj=%d, stats=%d, mem=%d, sparse=%f, rowscale=%d, ilut=%f)",
214 SetLevel_, SetBJ_, SetStats_, SetMem_, SetSparse_, SetRowScale_, SetILUT_);
215 // Set the parameters
216 eu->level = SetLevel_;
217 if(SetBJ_ != 0){
218 strcpy("bj", eu->algo_par);
219 }
220 if(SetSparse_ != 0.0){
221 eu->sparseTolA = SetSparse_;
222 }
223 if(SetRowScale_ != 0){
224 eu->isScaled = true;
225 }
226 if(SetILUT_ != 0.0){
227 eu->droptol = SetILUT_;
228 }
229 if(SetStats_ != 0 || SetMem_ != 0){
230 eu->logging = true;
231 Parser_dhInsert(parser_dh, "-eu_stats", "1");
232 }
233 // eu->A is the matrix as a void pointer, eu->m is local rows, eu->n is global rows
234 eu->A = (void*) A_.get();
235 eu->m = A_->NumMyRows();
236 eu->n = A_->NumGlobalRows();
237 Euclid_dhSetup(eu);
238 IsComputed_ = true;
239 NumCompute_ = NumCompute_ + 1;
240 ComputeTime_ = ComputeTime_ + Time_.ElapsedTime();
241 return 0;
242} //Compute()
243
244//==============================================================================
245int Ifpack_Euclid::ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const{
246 if(IsComputed() == false){
247 IFPACK_CHK_ERR(-1);
248 }
249 int NumVectors = X.NumVectors();
250 if(NumVectors != Y.NumVectors()){
251 IFPACK_CHK_ERR(-2);
252 }
253 Time_.ResetStartTime();
254 // Loop through the vectors
255 for(int vecNum = 0; vecNum < NumVectors; vecNum++){
256 CallEuclid(X[vecNum], Y[vecNum]);
257 }
258 if(SetStats_ != 0){
259 Euclid_dhPrintTestData(eu, stdout);
260 }
261 NumApplyInverse_ = NumApplyInverse_ + 1;
262 ApplyInverseTime_ = ApplyInverseTime_ + Time_.ElapsedTime();
263 return 0;
264} //ApplyInverse()
265
266//==============================================================================
267int Ifpack_Euclid::CallEuclid(double *x, double *y) const{
268 Euclid_dhApply(eu, x, y);
269 return 0;
270} //CallEuclid()
271
272//==============================================================================
273std::ostream& operator << (std::ostream& os, const Ifpack_Euclid& A){
274 if (!A.Comm().MyPID()) {
275 os << endl;
276 os << "================================================================================" << endl;
277 os << "Ifpack_Euclid: " << A.Label () << endl << endl;
278 os << "Using " << A.Comm().NumProc() << " processors." << endl;
279 os << "Global number of rows = " << A.Matrix().NumGlobalRows() << endl;
280 os << "Global number of nonzeros = " << A.Matrix().NumGlobalNonzeros() << endl;
281 os << "Condition number estimate = " << A.Condest() << endl;
282 os << endl;
283 os << "Phase # calls Total Time (s) Total MFlops MFlops/s" << endl;
284 os << "----- ------- -------------- ------------ --------" << endl;
285 os << "Initialize() " << std::setw(5) << A.NumInitialize()
286 << " " << std::setw(15) << A.InitializeTime()
287 << " 0.0 0.0" << endl;
288 os << "Compute() " << std::setw(5) << A.NumCompute()
289 << " " << std::setw(15) << A.ComputeTime()
290 << " " << std::setw(15) << 1.0e-6 * A.ComputeFlops();
291 if (A.ComputeTime() != 0.0)
292 os << " " << std::setw(15) << 1.0e-6 * A.ComputeFlops() / A.ComputeTime() << endl;
293 else
294 os << " " << std::setw(15) << 0.0 << endl;
295 os << "ApplyInverse() " << std::setw(5) << A.NumApplyInverse()
296 << " " << std::setw(15) << A.ApplyInverseTime()
297 << " " << std::setw(15) << 1.0e-6 * A.ApplyInverseFlops();
298 if (A.ApplyInverseTime() != 0.0)
299 os << " " << std::setw(15) << 1.0e-6 * A.ApplyInverseFlops() / A.ApplyInverseTime() << endl;
300 else
301 os << " " << std::setw(15) << 0.0 << endl;
302 os << "================================================================================" << endl;
303 os << endl;
304 }
305 return os;
306} // <<
307
308//==============================================================================
309double Ifpack_Euclid::Condest(const Ifpack_CondestType CT,
310 const int MaxIters,
311 const double Tol,
312 Epetra_RowMatrix* Matrix_in){
313 if (!IsComputed()) // cannot compute right now
314 return(-1.0);
315 return(Condest_);
316} //Condest() - not implemented
317
318#endif // HAVE_EUCLID && HAVE_MPI
int NumVectors() const