Amesos Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
example_AmesosFactory_Tridiag.cpp
Go to the documentation of this file.
1
2// @HEADER
3// ***********************************************************************
4//
5// Trilinos: An Object-Oriented Solver Framework
6// Copyright (2001) 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// This library is free software; you can redistribute it and/or modify
12// it under the terms of the GNU Lesser General Public License as
13// published by the Free Software Foundation; either version 2.1 of the
14// License, or (at your option) any later version.
15//
16// This library is distributed in the hope that it will be useful, but
17// WITHOUT ANY WARRANTY; without even the implied warranty of
18// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19// Lesser General Public License for more details.
20//
21// You should have received a copy of the GNU Lesser General Public
22// License along with this library; if not, write to the Free Software
23// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
24// USA
25// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
26//
27// ***********************************************************************
28// @HEADER
29
30#include "Amesos_ConfigDefs.h"
31#ifdef HAVE_MPI
32#include "mpi.h"
33#include "Epetra_MpiComm.h"
34#else
35#include "Epetra_SerialComm.h"
36#endif
37#include "Amesos_ConfigDefs.h"
38#include "Amesos.h"
39#include "Epetra_Map.h"
40#include "Epetra_CrsMatrix.h"
41#include "Epetra_Vector.h"
42#include "Epetra_LinearProblem.h"
43#include "Teuchos_ParameterList.hpp"
44
45// ==================== //
46// M A I N D R I V E R //
47// ==================== //
48//
49// This example will:
50// 1.- Create an tridiagonal matrix;
51// 2.- Call SymbolicFactorization();
52// 3.- Change the numerical values of the matrix;
53// 4.- Call NumericFactorization();
54// 5.- Set the entries of the RHS;
55// 6.- Call Solve().
56//
57// This example is intended to show the required data
58// for each phase. Phase (2) requires the matrix structure only.
59// Phase (4) requires the matrix structure (supposed unchanged
60// from phase (2)) and the matrix data. Phase (6) requires the
61// RHS and solution vector.
62//
63// This example can be run with any number of processors.
64//
65// Author: Marzio Sala, SNL 9214
66// Last modified: Apr-05.
67
68int main(int argc, char *argv[])
69{
70
71#ifdef HAVE_MPI
72 MPI_Init(&argc, &argv);
73 Epetra_MpiComm Comm(MPI_COMM_WORLD);
74#else
75 Epetra_SerialComm Comm;
76#endif
77
78 int NumGlobalElements = 100; // global dimension of the problem.
79
80 // ======================================================= //
81 // B E G I N N I N G O F M A T R I X C R E A T I O N //
82 // ======================================================= //
83
84 // Construct a Map that puts approximatively the same number of
85 // equations on each processor. `0' is the index base (that is,
86 // numbering starts from 0.
87 Epetra_Map Map(NumGlobalElements, 0, Comm);
88
89 // Create an empty EpetraCrsMatrix
90 Epetra_CrsMatrix A(Copy, Map, 0);
91
92 // Create the structure of the matrix (tridiagonal)
93 int NumMyElements = Map.NumMyElements();
94
95 // Add rows one-at-a-time
96 // Need some vectors to help
97
98 double Values[3];
99 // Right now, we put zeros only in the matrix.
100 Values[0] = 0.0;
101 Values[1] = 0.0;
102 Values[2] = 0.0;
103 int Indices[3];
104 int NumEntries;
106 int* MyGlobalElements = Map.MyGlobalElements();
107
108 // At this point we simply set the nonzero structure of A.
109 // Actual values will be inserted later (now all zeros)
110 for (int i = 0; i < NumMyElements; i++)
111 {
112 if (MyGlobalElements[i] == 0)
113 {
114 Indices[0] = 0;
115 Indices[1] = 1;
116 NumEntries = 2;
117 }
118 else if (MyGlobalElements[i] == NumGlobalElements-1)
119 {
120 Indices[0] = NumGlobalElements-1;
121 Indices[1] = NumGlobalElements-2;
122 NumEntries = 2;
123 }
124 else
125 {
126 Indices[0] = MyGlobalElements[i]-1;
127 Indices[1] = MyGlobalElements[i];
128 Indices[2] = MyGlobalElements[i]+1;
129 NumEntries = 3;
130 }
131
132 AMESOS_CHK_ERR(A.InsertGlobalValues(MyGlobalElements[i],
133 NumEntries, Values, Indices));
134 }
135
136 // Finish up.
137 A.FillComplete();
138
139 // =========================================== //
140 // E N D O F M A T R I X C R E A T I O N //
141 // =========================================== //
142
143 // Now the matrix STRUCTURE is set. We cannot add
144 // new nonzero elements, but we can still change the
145 // numerical values of all inserted elements (as we will
146 // do later).
147
148 // ===================================================== //
149 // B E G I N N I N G O F T H E AM E S O S P A R T //
150 // ===================================================== //
151
152 // For comments on the commands in this section, please
153 // see file example_AmesosFactory.cpp.
154
155 Epetra_LinearProblem Problem;
156
157 Problem.SetOperator(&A);
158
159 // Initializes Amesos solver. Here we solve with Amesos_Klu.
160
161 Amesos_BaseSolver * Solver;
162 Amesos Factory;
163 Solver = Factory.Create("Amesos_Klu", Problem);
164
165 // Factory.Create() returns 0 if the requested solver
166 // is not available
167
168 if (Solver == 0) {
169 std::cerr << "Selected solver is not available" << std::endl;
170 // return ok not to break the test harness
171#ifdef HAVE_MPI
172 MPI_Finalize();
173#endif
174 exit(EXIT_SUCCESS);
175 }
176
177 // At this point we can perform the numeric factorization.
178 // Note that the matrix contains 0's only.
179
180 Solver->SymbolicFactorization();
181
182 // Now, we repopulate the matrix with entries corresponding
183 // to a 1D Laplacian. LHS and RHS are still untouched.
184
185 for (int i = 0; i < NumMyElements; i++)
186 {
187 if (MyGlobalElements[i] == 0)
188 {
189 Indices[0] = 0;
190 Indices[1] = 1;
191 Values[0] = 2.0;
192 Values[1] = -1.0;
193 NumEntries = 2;
194 }
195 else if (MyGlobalElements[i] == NumGlobalElements-1)
196 {
197 Indices[0] = NumGlobalElements - 1;
198 Indices[1] = NumGlobalElements - 2;
199 Values[0] = 2.0;
200 Values[1] = -1.0;
201 NumEntries = 2;
202 }
203 else
204 {
205 Indices[0] = MyGlobalElements[i] - 1;
206 Indices[1] = MyGlobalElements[i];
207 Indices[2] = MyGlobalElements[i] + 1;
208 Values[0] = -1.0;
209 Values[1] = 2.0;
210 Values[2] = -1.0;
211 NumEntries = 3;
212 }
213
214 AMESOS_CHK_ERR(A.ReplaceGlobalValues(MyGlobalElements[i],
215 NumEntries, Values, Indices));
216 }
217
218 // ... and we can compute the numeric factorization.
219 Solver->NumericFactorization();
220
221 // Finally, we set up the LHS and the RHS vector (Random().
222 Epetra_Vector b(Map);
223 b.Random();
224 Epetra_Vector x(Map);
225 x.PutScalar(0.0);
226
227 Problem.SetLHS(&x);
228 Problem.SetRHS(&b);
229
230 Solver->Solve();
231
232 // Print out the timing information and get it from the solver
233 Solver->PrintTiming();
234
235 Teuchos::ParameterList TimingsList;
236 Solver->GetTiming( TimingsList );
237
238 // You can find out how much time was spent in ...
239 double sfact_time, nfact_time, solve_time;
240 double mtx_conv_time, mtx_redist_time, vec_redist_time;
241
242 // 1) The symbolic factorization
243 // (parameter doesn't always exist)
244 sfact_time = TimingsList.get( "Total symbolic factorization time", 0.0 );
245
246 // 2) The numeric factorization
247 // (always exists if NumericFactorization() is called)
248 nfact_time = Teuchos::getParameter<double>( TimingsList, "Total numeric factorization time" );
249
250 // 3) Solving the linear system
251 // (always exists if Solve() is called)
252 solve_time = Teuchos::getParameter<double>( TimingsList, "Total solve time" );
253
254 // 4) Converting the matrix to the accepted format for the solver
255 // (always exists if SymbolicFactorization() is called)
256 mtx_conv_time = Teuchos::getParameter<double>( TimingsList, "Total solve time" );
257
258 // 5) Redistributing the matrix for each solve to the accepted format for the solver
259 mtx_redist_time = TimingsList.get( "Total matrix redistribution time", 0.0 );
260
261 // 6) Redistributing the vector for each solve to the accepted format for the solver
262 vec_redist_time = TimingsList.get( "Total vector redistribution time", 0.0 );
263
264 std::cout << " sfact_time " << sfact_time
265 << " nfact_time " << nfact_time
266 << " solve_time " << solve_time
267 << " mtx_conv_time " << mtx_conv_time
268 << " mtx_redist_time " << mtx_redist_time
269 << " vec_redist_time " << vec_redist_time
270 << std::endl;
271
272 // =========================================== //
273 // E N D O F T H E A M E S O S P A R T //
274 // =========================================== //
275
276 // ================== //
277 // compute ||Ax - b|| //
278 // ================== //
279
280 double residual;
281
282 Epetra_Vector Ax(b.Map());
283 A.Multiply(false, x, Ax);
284 Ax.Update(1.0, b, -1.0);
285 Ax.Norm2(&residual);
286
287 if (!Comm.MyPID())
288 std::cout << "After AMESOS solution, ||b-Ax||_2 = " << residual << std::endl;
289
290 // delete Solver. Do this before MPI_Finalize()
291 // as MPI calls can occur in the destructor.
292 delete Solver;
293
294 if (residual > 1e-5)
295 return(EXIT_FAILURE);
296
297#ifdef HAVE_MPI
298 MPI_Finalize();
299#endif
300 return(EXIT_SUCCESS);
301
302} // end of main()
#define AMESOS_CHK_ERR(a)
Amesos_BaseSolver: A pure virtual class for direct solution of real-valued double-precision operators...
virtual void PrintTiming() const =0
Prints timing information about the current solver.
virtual int Solve()=0
Solves A X = B (or AT x = B)
virtual void GetTiming(Teuchos::ParameterList &TimingParameterList) const
Extracts timing information from the current solver and places it in the parameter list....
virtual int NumericFactorization()=0
Performs NumericFactorization on the matrix A.
virtual int SymbolicFactorization()=0
Performs SymbolicFactorization on the matrix A.
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
int main(int argc, char *argv[])