← Back to branch summary

~ubuntu-branches/ubuntu/wily/dolfin/wily-proposed

« back to all changes in this revision

Viewing changes to dolfin/la/EpetraLUSolver.cpp

  • Committer: Package Import Robot
  • Author(s): Johannes Ring
  • Date: 2015-03-17 07:57:11 UTC
  • mfrom: (1.1.18) (19.1.24 experimental)
  • Revision ID: package-import@ubuntu.com-20150317075711-1v207zbty9qmygow
Tags: 1.5.0-1
http://bugs.debian.org/780359
* New upstream release (closes: #780359).
* debian/control:
  - Bump Standards-Version to 3.9.6 (no changes needed).
  - Bump X-Python-Version to >= 2.7.
  - Update package names for new SONAME 1.5 (libdolfin1.4 ->
    libdolfin1.5, libdolfin1.4-dbg -> libdolfin1.5-dbg and
    libdolfin1.4-dev -> libdolfin1.5-dev).
  - Bump minimum required version for python-instant, python-ufl and
    python-ffc to 1.5.0.
  - Add python-sympy and python-six to Depends for binary package
    python-dolfin.
  - Add dh-python to Build-Depends.
  - Remove libcgal-dev from {Build-}Depends.
* Remove CSGCGALMeshGenerator3D-oom.patch since CGAL is no longer used
  by DOLFIN.
* Move debian/libdolfin1.4.install -> debian/libdolfin1.5.install.
* debian/rules: No longer any non DFSG-free stuff to remove, so update
  get-orig-source target (update debian/watch accordingly).
* Update debian/copyright and debian/copyright_hints.

Show diffs side-by-side

added added

removed removed

Lines of Context:
dolfin/la/EpetraLUSolver.cpp
1
 
// Copyright (C) 2008-2010 Kent-Andre Mardal and Garth N. Wells
2
 
//
3
 
// This file is part of DOLFIN.
4
 
//
5
 
// DOLFIN is free software: you can redistribute it and/or modify
6
 
// it under the terms of the GNU Lesser General Public License as published by
7
 
// the Free Software Foundation, either version 3 of the License, or
8
 
// (at your option) any later version.
9
 
//
10
 
// DOLFIN is distributed in the hope that it will be useful,
11
 
// but WITHOUT ANY WARRANTY; without even the implied warranty of
12
 
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13
 
// GNU Lesser General Public License for more details.
14
 
//
15
 
// You should have received a copy of the GNU Lesser General Public License
16
 
// along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
17
 
//
18
 
// Modified by Anders Logg 2011-2012
19
 
//
20
 
// First added:  2008
21
 
// Last changed: 2012-08-22
22
 
 
23
 
#ifdef HAS_TRILINOS
24
 
 
25
 
// Included here to avoid a C++ problem with some MPI implementations
26
 
#include <dolfin/common/MPI.h>
27
 
 
28
 
#include <Amesos.h>
29
 
#include <Amesos_BaseSolver.h>
30
 
#include <Amesos_ConfigDefs.h>
31
 
#include <Amesos_Klu.h>
32
 
#include <Epetra_FECrsMatrix.h>
33
 
#include <Epetra_FEVector.h>
34
 
#include <Epetra_LinearProblem.h>
35
 
#include <Epetra_MultiVector.h>
36
 
 
37
 
#include <dolfin/common/MPI.h>
38
 
#include <dolfin/common/NoDeleter.h>
39
 
#include <dolfin/common/Timer.h>
40
 
#include "GenericLinearOperator.h"
41
 
#include "GenericVector.h"
42
 
#include "EpetraMatrix.h"
43
 
#include "EpetraVector.h"
44
 
#include "LUSolver.h"
45
 
#include "EpetraLUSolver.h"
46
 
 
47
 
using namespace dolfin;
48
 
 
49
 
//-----------------------------------------------------------------------------
50
 
std::vector<std::pair<std::string, std::string> >
51
 
EpetraLUSolver::methods()
52
 
{
53
 
  static std::vector<std::pair<std::string, std::string> > m;
54
 
 
55
 
  m.push_back(std::make_pair("default", "default LU solver"));
56
 
  m.push_back(std::make_pair("umfpack", "UMFPACK (Unsymmetric MultiFrontal sparse LU factorization)"));
57
 
  m.push_back(std::make_pair("mumps", "MUMPS (MUltifrontal Massively Parallel Sparse direct Solver)"));
58
 
  m.push_back(std::make_pair("klu", "Trilinos KLU"));
59
 
 
60
 
  return m;
61
 
}
62
 
//-----------------------------------------------------------------------------
63
 
std::string EpetraLUSolver::choose_method(std::string method) const
64
 
{
65
 
  Amesos factory;
66
 
  if (method == "default")
67
 
  {
68
 
    if (factory.Query("Amesos_Mumps"))
69
 
      method = "Amesos_Mumps";
70
 
    else if (factory.Query("Amesos_Umfpack"))
71
 
      method = "Amesos_Umfpack";
72
 
    else if (factory.Query("Amesos_Klu"))
73
 
      method = "Amesos_Klu";
74
 
    else
75
 
    {
76
 
      dolfin_error("EpetraLUSolver.cpp",
77
 
                   "choose default Epetra LU solver",
78
 
                   "No methods available");
79
 
    }
80
 
  }
81
 
  else if (method == "umfpack")
82
 
    method = "Amesos_Umfpack";
83
 
  else if (method == "mumps")
84
 
    method = "Amesos_mumps";
85
 
  else if (method == "klu")
86
 
    method = "Amesos_Klu";
87
 
  else
88
 
  {
89
 
    dolfin_error("EpetraLUSolver.cpp",
90
 
                 "solve linear system",
91
 
                 "Unknown LU solver method \"%s\". "
92
 
                 "Use list_lu_solver_methods() to list available methods",
93
 
                 method.c_str());
94
 
  }
95
 
 
96
 
  return method;
97
 
}
98
 
//-----------------------------------------------------------------------------
99
 
Parameters EpetraLUSolver::default_parameters()
100
 
{
101
 
  Parameters p(LUSolver::default_parameters());
102
 
  p.rename("epetra_lu_solver");
103
 
  return p;
104
 
}
105
 
//-----------------------------------------------------------------------------
106
 
EpetraLUSolver::EpetraLUSolver(std::string method)
107
 
  : symbolic_factorized(false),
108
 
    numeric_factorized(false),
109
 
    linear_problem(new Epetra_LinearProblem)
110
 
{
111
 
  // Set default parameters
112
 
  parameters = default_parameters();
113
 
 
114
 
  // Choose method
115
 
  _method = choose_method(method);
116
 
 
117
 
  // Initialize solver
118
 
  Amesos factory;
119
 
  solver.reset(factory.Create(_method, *linear_problem));
120
 
 
121
 
  // Check that solver was initialized correctly
122
 
  if (!solver)
123
 
  {
124
 
    dolfin_error("EpetraLUSolver.cpp",
125
 
                 "create Epetra LU solver",
126
 
                 "Epetra was not able to create linear solver \"%s\"",
127
 
                 _method.c_str());
128
 
  }
129
 
}
130
 
//-----------------------------------------------------------------------------
131
 
EpetraLUSolver::EpetraLUSolver(std::shared_ptr<const GenericLinearOperator> A,
132
 
                               std::string method)
133
 
  : symbolic_factorized(false),
134
 
    numeric_factorized(false),
135
 
    linear_problem(new Epetra_LinearProblem)
136
 
{
137
 
  // Set default parameters
138
 
  parameters = default_parameters();
139
 
 
140
 
  // Set operator
141
 
  _matA = as_type<const EpetraMatrix>(require_matrix(A));
142
 
  dolfin_assert(_matA);
143
 
 
144
 
  // Choose method
145
 
  _method = choose_method(method);
146
 
 
147
 
  // Initialize solver
148
 
  Amesos factory;
149
 
  solver.reset(factory.Create(_method, *linear_problem));
150
 
 
151
 
  // Check that solver was initialized correctly
152
 
  if (!solver)
153
 
  {
154
 
    dolfin_error("EpetraLUSolver.cpp",
155
 
                 "create Epetra LU solver",
156
 
                 "Epetra was not able to create linear solver \"%s\"",
157
 
                 _method.c_str());
158
 
  }
159
 
}
160
 
//-----------------------------------------------------------------------------
161
 
EpetraLUSolver::~EpetraLUSolver()
162
 
{
163
 
  // Do nothing
164
 
}
165
 
//-----------------------------------------------------------------------------
166
 
void
167
 
EpetraLUSolver::set_operator(std::shared_ptr<const GenericLinearOperator> A)
168
 
{
169
 
  dolfin_assert(linear_problem);
170
 
 
171
 
  _matA = as_type<const EpetraMatrix>(require_matrix(A));
172
 
  dolfin_assert(_matA);
173
 
  linear_problem->SetOperator(_matA->mat().get());
174
 
 
175
 
  symbolic_factorized = false;
176
 
  numeric_factorized  = false;
177
 
}
178
 
//-----------------------------------------------------------------------------
179
 
const GenericLinearOperator& EpetraLUSolver::get_operator() const
180
 
{
181
 
  if (!_matA)
182
 
  {
183
 
    dolfin_error("EpetraLUSolver.cpp",
184
 
                 "access operator for Epetra LU solver",
185
 
                 "Operator has not been set");
186
 
  }
187
 
  return *_matA;
188
 
}
189
 
//-----------------------------------------------------------------------------
190
 
std::size_t EpetraLUSolver::solve(GenericVector& x, const GenericVector& b)
191
 
{
192
 
  Timer timer("Epetra LU solver");
193
 
 
194
 
  dolfin_assert(linear_problem);
195
 
  dolfin_assert(solver);
196
 
 
197
 
  // Write a message
198
 
  if (parameters["report"] && dolfin::MPI::rank(MPI_COMM_WORLD) == 0)
199
 
  {
200
 
    if (solver->UseTranspose())
201
 
    {
202
 
      info("Solving linear system transposed of size %d x %d using Epetra LU solver (%s).",
203
 
           _matA->size(0), _matA->size(1), _method.c_str());
204
 
    }
205
 
    else
206
 
    {
207
 
      info("Solving linear system of size %d x %d using Epetra LU solver (%s).",
208
 
           _matA->size(0), _matA->size(1), _method.c_str());
209
 
    }
210
 
  }
211
 
 
212
 
  // Downcast vector
213
 
  EpetraVector& _x = as_type<EpetraVector>(x);
214
 
  const EpetraVector& _b = as_type<const EpetraVector>(b);
215
 
 
216
 
  // Get operator matrix
217
 
  const Epetra_RowMatrix* A = linear_problem->GetMatrix();
218
 
  if (!A)
219
 
  {
220
 
    dolfin_error("EpetraLUSolver.cpp",
221
 
                 "solve linear system using Epetra LU solver",
222
 
                 "Operator has not been set");
223
 
  }
224
 
 
225
 
  const std::size_t N = A->NumGlobalCols64();
226
 
  if (N != b.size())
227
 
  {
228
 
    dolfin_error("EpetraLUSolver.cpp",
229
 
                 "solve linear system using Epetra LU solver",
230
 
                 "Non-matching dimensions for linear system");
231
 
  }
232
 
 
233
 
  // Initialize solution vector
234
 
  if (x.empty())
235
 
  {
236
 
    _matA->init_vector(x, 1);
237
 
    x.zero();
238
 
  }
239
 
 
240
 
  // Set LHS and RHS vectors
241
 
  linear_problem->SetRHS(_b.vec().get());
242
 
  linear_problem->SetLHS(_x.vec().get());
243
 
 
244
 
  // Get some parameters
245
 
  const bool reuse_fact   = parameters["reuse_factorization"];
246
 
  const bool same_pattern = parameters["same_nonzero_pattern"];
247
 
 
248
 
  // Perform symbolic factorization
249
 
  if ( (reuse_fact || same_pattern) && !symbolic_factorized )
250
 
  {
251
 
    AMESOS_CHK_ERR(solver->SymbolicFactorization());
252
 
    symbolic_factorized = true;
253
 
  }
254
 
  else if (!reuse_fact && !same_pattern)
255
 
  {
256
 
    AMESOS_CHK_ERR(solver->SymbolicFactorization());
257
 
    symbolic_factorized = true;
258
 
  }
259
 
 
260
 
  // Perform numeric factorization
261
 
  if (reuse_fact && !numeric_factorized)
262
 
  {
263
 
    AMESOS_CHK_ERR(solver->NumericFactorization());
264
 
    numeric_factorized = true;
265
 
  }
266
 
  else if (!reuse_fact)
267
 
  {
268
 
    AMESOS_CHK_ERR(solver->NumericFactorization());
269
 
    numeric_factorized = true;
270
 
  }
271
 
 
272
 
  // Solve
273
 
  AMESOS_CHK_ERR(solver->Solve());
274
 
 
275
 
  // Update ghost values
276
 
  _x.update_ghost_values();
277
 
 
278
 
  return 1;
279
 
}
280
 
//-----------------------------------------------------------------------------
281
 
std::size_t EpetraLUSolver::solve(const GenericLinearOperator& A,
282
 
                                  GenericVector& x,
283
 
                                  const GenericVector& b)
284
 
{
285
 
  return solve(as_type<const EpetraMatrix>(require_matrix(A)),
286
 
               as_type<EpetraVector>(x),
287
 
               as_type<const EpetraVector>(b));
288
 
}
289
 
//-----------------------------------------------------------------------------
290
 
std::size_t EpetraLUSolver::solve(const EpetraMatrix& A, EpetraVector& x,
291
 
                                  const EpetraVector& b)
292
 
{
293
 
  std::shared_ptr<const EpetraMatrix> Atmp(&A, NoDeleter());
294
 
  set_operator(Atmp);
295
 
  return solve(x, b);
296
 
}
297
 
//-----------------------------------------------------------------------------
298
 
std::size_t EpetraLUSolver::solve_transpose(GenericVector& x,
299
 
                                            const GenericVector& b)
300
 
{
301
 
  dolfin_assert(solver);
302
 
  solver->SetUseTranspose(true);
303
 
  std::size_t out = solve(x, b);
304
 
  solver->SetUseTranspose(false);
305
 
  return out;
306
 
}
307
 
//-----------------------------------------------------------------------------
308
 
std::size_t EpetraLUSolver::solve_transpose(const GenericLinearOperator& A,
309
 
                                            GenericVector& x,
310
 
                                            const GenericVector& b)
311
 
{
312
 
  dolfin_assert(solver);
313
 
  solver->SetUseTranspose(true);
314
 
  std::size_t out = solve(A, x, b);
315
 
  solver->SetUseTranspose(false);
316
 
  return out;
317
 
}
318
 
//-----------------------------------------------------------------------------
319
 
std::size_t EpetraLUSolver::solve_transpose(const EpetraMatrix& A,
320
 
                                            EpetraVector& x,
321
 
                                            const EpetraVector& b)
322
 
{
323
 
  dolfin_assert(solver);
324
 
  solver->SetUseTranspose(true);
325
 
  std::size_t out = solve(A, x, b);
326
 
  solver->SetUseTranspose(false);
327
 
  return out;
328
 
}
329
 
//-----------------------------------------------------------------------------
330
 
std::string EpetraLUSolver::str(bool verbose) const
331
 
{
332
 
  dolfin_not_implemented();
333
 
  return std::string();
334
 
}
335
 
//-----------------------------------------------------------------------------
336
 
 
337
 
#endif