Amesos Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
RunParaklete.cpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Amesos: An Interface to Direct Solvers
5// Copyright (2004) 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// This library is free software; you can redistribute it and/or modify
11// it under the terms of the GNU Lesser General Public License as
12// published by the Free Software Foundation; either version 2.1 of the
13// License, or (at your option) any later version.
14//
15// This library is distributed in the hope that it will be useful, but
16// WITHOUT ANY WARRANTY; without even the implied warranty of
17// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18// Lesser General Public License for more details.
19//
20// You should have received a copy of the GNU Lesser General Public
21// License along with this library; if not, write to the Free Software
22// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23// USA
24// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
25//
26// ***********************************************************************
27// @HEADER
28
29#include "Amesos_ConfigDefs.h"
30// This example needs Galeri to generate the linear system.
31#ifdef HAVE_MPI
32#include "mpi.h"
33#include "Epetra_MpiComm.h"
34#else
35#include "Epetra_SerialComm.h"
36#endif
37#include "Epetra_Vector.h"
38#include "Epetra_Time.h"
39#include "Epetra_RowMatrix.h"
40#include "Epetra_CrsMatrix.h"
41#include "Amesos.h"
42#include "Amesos_BaseSolver.h"
43#include "Teuchos_ParameterList.hpp"
44#include "Galeri_Maps.h"
45#include "Galeri_CrsMatrices.h"
46
47using namespace Teuchos;
48using namespace Galeri;
49
50int MyCreateCrsMatrix( const char *in_filename, const Epetra_Comm &Comm,
51 Epetra_Map *& readMap,
52 const bool transpose, const bool distribute,
53 bool& symmetric, Epetra_CrsMatrix *& Matrix ) {
54
55 Epetra_CrsMatrix * readA = 0;
56 Epetra_Vector * readx = 0;
57 Epetra_Vector * readb = 0;
58 Epetra_Vector * readxexact = 0;
59
60 //
61 // This hack allows TestOptions to be run from either the test/TestOptions/ directory or from
62 // the test/ directory (as it is in nightly testing and in make "run-tests")
63 //
64 FILE *in_file = fopen( in_filename, "r");
65
66 char *filename;
67 if (in_file == NULL )
68 filename = &in_filename[1] ; // Strip off ithe "." from
69 // "../" and try again
70 else {
71 filename = in_filename ;
72 fclose( in_file );
73 }
74
75 symmetric = false ;
76 std::string FileName = filename ;
77
78 int FN_Size = FileName.size() ;
79 std::string LastFiveBytes = FileName.substr( EPETRA_MAX(0,FN_Size-5), FN_Size );
80 std::string LastFourBytes = FileName.substr( EPETRA_MAX(0,FN_Size-4), FN_Size );
81
82 if ( LastFiveBytes == ".triU" ) {
83 // Call routine to read in unsymmetric Triplet matrix
84 EPETRA_CHK_ERR( Trilinos_Util_ReadTriples2Epetra( filename, false, Comm, readMap, readA, readx,
85 readb, readxexact) );
86 symmetric = false;
87 } else {
88 if ( LastFiveBytes == ".triS" ) {
89 // Call routine to read in symmetric Triplet matrix
90 EPETRA_CHK_ERR( Trilinos_Util_ReadTriples2Epetra( filename, true, Comm, readMap, readA, readx,
91 readb, readxexact) );
92 symmetric = true;
93 } else {
94 if ( LastFourBytes == ".mtx" ) {
95 EPETRA_CHK_ERR( Trilinos_Util_ReadMatrixMarket2Epetra( filename, Comm, readMap,
96 readA, readx, readb, readxexact) );
97 FILE* in_file = fopen( filename, "r");
98 assert (in_file != NULL) ; // Checked in Trilinos_Util_CountMatrixMarket()
99 const int BUFSIZE = 800 ;
100 char buffer[BUFSIZE] ;
101 fgets( buffer, BUFSIZE, in_file ) ; // Pick symmetry info off of this string
102 std::string headerline1 = buffer;
103#ifdef TFLOP
104 if ( headerline1.find("symmetric") < BUFSIZE ) symmetric = true;
105#else
106 if ( headerline1.find("symmetric") != std::string::npos) symmetric = true;
107
108#endif
109 fclose(in_file);
110
111 } else {
112 // Call routine to read in HB problem
113 Trilinos_Util_ReadHb2Epetra( filename, Comm, readMap, readA, readx,
114 readb, readxexact) ;
115 if ( LastFourBytes == ".rsa" ) symmetric = true ;
116 }
117 }
118 }
119
120
121 if ( readb ) delete readb;
122 if ( readx ) delete readx;
123 if ( readxexact ) delete readxexact;
124
125 Epetra_CrsMatrix *serialA ;
126 Epetra_CrsMatrix *transposeA;
127
128 if ( transpose ) {
129 transposeA = new Epetra_CrsMatrix( Copy, *readMap, 0 );
130 assert( CrsMatrixTranspose( readA, transposeA ) == 0 );
131 serialA = transposeA ;
132 delete readA;
133 readA = 0 ;
134 } else {
135 serialA = readA ;
136 }
137
138 assert( (void *) &serialA->Graph() ) ;
139 assert( (void *) &serialA->RowMap() ) ;
140 assert( serialA->RowMap().SameAs(*readMap) ) ;
141
142 if ( distribute ) {
143 // Create uniform distributed map
144 Epetra_Map DistMap(readMap->NumGlobalElements(), 0, Comm);
145
146 // Create Exporter to distribute read-in matrix and vectors
147 Epetra_Export exporter( *readMap, DistMap );
148
149 Epetra_CrsMatrix *Amat = new Epetra_CrsMatrix( Copy, DistMap, 0 );
150 Amat->Export(*serialA, exporter, Add);
151 assert(Amat->FillComplete()==0);
152
153 Matrix = Amat;
154 //
155 // Make sure that deleting Amat->RowMap() will delete map
156 //
157 // Bug: We can't manage to delete map his way anyway,
158 // and this fails on tranposes, so for now I just accept
159 // the memory loss.
160 // assert( &(Amat->RowMap()) == map ) ;
161 delete readMap;
162 readMap = 0 ;
163 delete serialA;
164 } else {
165
166 Matrix = serialA;
167 }
168
169
170 return 0;
171}
172
173
174
175
176
177
178// ===================== //
179// M A I N D R I V E R //
180// ===================== //
181//
182// This example compares all the available Amesos solvers
183// for the solution of the same linear system.
184//
185// The example can be run in serial and in parallel.
186//
187// Author: Marzio Sala, SNL 9214
188// Last modified: Oct-05.
189
190int main(int argc, char *argv[])
191{
192#ifdef HAVE_MPI
193 MPI_Init(&argc, &argv);
194 Epetra_MpiComm Comm(MPI_COMM_WORLD);
195#else
196 Epetra_SerialComm Comm;
197#endif
198
199 bool verbose = (Comm.MyPID() == 0);
200 double TotalResidual = 0.0;
201
202 // Create the Map, defined as a grid, of size nx x ny x nz,
203 // subdivided into mx x my x mz cubes, each assigned to a
204 // different processor.
205
206#if 0
207 ParameterList GaleriList;
208 GaleriList.set("nx", 4);
209 GaleriList.set("ny", 4);
210 GaleriList.set("nz", 4 * Comm.NumProc());
211 GaleriList.set("mx", 1);
212 GaleriList.set("my", 1);
213 GaleriList.set("mz", Comm.NumProc());
214 Epetra_Map* Map = CreateMap("Cartesian3D", Comm, GaleriList);
215
216 // Create a matrix, in this case corresponding to a 3D Laplacian
217 // discretized using a classical 7-point stencil. Please refer to
218 // the Galeri documentation for an overview of available matrices.
219 //
220 // NOTE: matrix must be symmetric if DSCPACK is used.
221
222 Epetra_CrsMatrix* Matrix = CreateCrsMatrix("Laplace3D", Map, GaleriList);
223#else
224 bool transpose = false ;
225 bool distribute = false ;
226 bool symmetric ;
227 Epetra_CrsMatrix *Matrix = 0 ;
228 Epetra_Map *Map = 0 ;
229 CreateCrsMatrix( "ibm.triU", Comm, Map, transpose, distribute, &symmetric, Matrix ) ;
230
231
232
233
234#endif
235
236 // build vectors, in this case with 1 vector
237 Epetra_MultiVector LHS(*Map, 1);
238 Epetra_MultiVector RHS(*Map, 1);
239
240 // create a linear problem object
241 Epetra_LinearProblem Problem(Matrix, &LHS, &RHS);
242
243 // use this list to set up parameters, now it is required
244 // to use all the available processes (if supported by the
245 // underlying solver). Uncomment the following two lines
246 // to let Amesos print out some timing and status information.
247 ParameterList List;
248 List.set("PrintTiming",true);
249 List.set("PrintStatus",true);
250 List.set("MaxProcs",Comm.NumProc());
251
252 std::vector<std::string> SolverType;
253 SolverType.push_back("Amesos_Paraklete");
254
255 Epetra_Time Time(Comm);
256
257 // this is the Amesos factory object that will create
258 // a specific Amesos solver.
259 Amesos Factory;
260
261 // Cycle over all solvers.
262 // Only installed solvers will be tested.
263 for (unsigned int i = 0 ; i < SolverType.size() ; ++i)
264 {
265 // Check whether the solver is available or not
266 if (Factory.Query(SolverType[i]))
267 {
268 // 1.- set exact solution (constant vector)
269 LHS.PutScalar(1.0);
270
271 // 2.- create corresponding rhs
272 Matrix->Multiply(false, LHS, RHS);
273
274 // 3.- randomize solution vector
275 LHS.Random();
276
277 // 4.- create the amesos solver object
278 Amesos_BaseSolver* Solver = Factory.Create(SolverType[i], Problem);
279 assert (Solver != 0);
280
281 Solver->SetParameters(List);
282
283 // 5.- factorize and solve
284
285 Time.ResetStartTime();
287
288 // 7.- delete the object
289 delete Solver;
290
291 }
292 }
293
294 delete Matrix;
295 delete Map;
296
297 if (TotalResidual > 1e-9)
298 exit(EXIT_FAILURE);
299
300#ifdef HAVE_MPI
301 MPI_Finalize();
302#endif
303
304 return(EXIT_SUCCESS);
305} // end of main()
static bool verbose
Definition: Amesos.cpp:67
#define AMESOS_CHK_ERR(a)
int CrsMatrixTranspose(Epetra_CrsMatrix *In, Epetra_CrsMatrix *Out)
int main(int argc, char *argv[])
int MyCreateCrsMatrix(const char *in_filename, const Epetra_Comm &Comm, Epetra_Map *&readMap, const bool transpose, const bool distribute, bool &symmetric, Epetra_CrsMatrix *&Matrix)
#define EPETRA_CHK_ERR(xxx)
int CreateCrsMatrix(const char *in_filename, const Epetra_Comm &Comm, Epetra_Map *&readMap, const bool transpose, const bool distribute, bool &symmetric, Epetra_CrsMatrix *&Matrix)
Definition: TestOptions.cpp:81
Amesos_BaseSolver: A pure virtual class for direct solution of real-valued double-precision operators...
virtual int SymbolicFactorization()=0
Performs SymbolicFactorization on the matrix A.
virtual int SetParameters(Teuchos::ParameterList &ParameterList)=0
Updates internal variables.
Factory for binding a third party direct solver to an Epetra_LinearProblem.
Definition: Amesos.h:44
Amesos_BaseSolver * Create(const char *ClassType, const Epetra_LinearProblem &LinearProblem)
Amesos Create method.
Definition: Amesos.cpp:69
bool Query(const char *ClassType)
Queries whether a given interface is available or not.
Definition: Amesos.cpp:193