dolfin team mailing list archive
-
dolfin team
-
Mailing list archive
-
Message #21064
Re: Some patches for dolfin (resend plus new)
As promised elsewhere, here is a patch which fixes parallel matrix
multiplication (#12). I've also rebased the old ones on top of the current
bzr version, and a new one. The description of patches up to 0009 is quoted
below, the new one is:
10) Replace all checks for isinstance(x, (float,int)) with numpy.isscalar(x)
in the python/swig layer. The reason is that numpy.float64 for example
should be supported.
-j.
On 20 January 2011 11:29, Joachim Berdal Haga <jobh@xxxxxxxxx> wrote:
> I've attached some patches for dolfin. Please apply (or not). Short
> description follows, and a question...
>
> 1) Allow tensor Constants in python (initialised from e.g. a numpy matrix)
>
> 2) Define a split method on FunctionSpace in python. This seems natural
> since Function has one... but not important.
>
> 3) Define the methods "near" and "between" in both compiled and python
> Expressions (etc). So you can write
> def boundary(self):
> return near(x[0], 0)
> or
> boundary = compile_subdomains("near(x[0], 0)")
>
> I think "near" is good and gets rid of many of those shouting
> DOLFIN_EPSs. "between" is perhaps more questionable... I leave it to you.
>
> 4) This matches up with a viper patch sent to fenics-viper@lp, and
> simplifies plot() a bit. It also does less projection of those types viper
> can actually handle --- which means that for example DG plots are slightly
> different (but maybe not better?).
>
> 5) I find it easier to understand DIVERGED_INDEFINITE_MAT than -8 as an
> error message.
>
> 6) PETSc development version requires some changes
>
> 7) This one is really unimportant and you may not agree anyway.
>
> 8) This one IS important. Allows assembly into rectangular matrices with
> epetra. This used to crash, and tracing it down was a bitch due to epetra
> throwing bare-int exceptions (which are rethrown by swig).
>
> 9) A utility function on DirichletBC for symmetric modification of block
> matrices. Also implements EpetraMatrix::set_row. These can probably be made
> much more efficient by someone who knows the la backends better (I threw in
> a "this is inefficient" log msg there to make this point), but it seems to
> work in superficial testing.
>
> Now the question: I seem to hit dolfin asserts a lot when I try things, and
> that's annoying since they call abort() and thus kill (i)python. I see that
> they used to be implemented as runtime exceptions, which plays much better
> with python, but that was changed. Was there a good reason for that change,
> or can we bring back exceptions? (the change itself is trivial and I can
> send a patch, but there was probably a reason...)
>
> --
> -j.
>
From 8dacd49b9e73670b201a80ea205acd27b002bb87 Mon Sep 17 00:00:00 2001
From: Joachim B Haga <jobh@xxxxxxxxx>
Date: Fri, 28 Jan 2011 14:19:50 +0100
Subject: [PATCH 9/9] Introduce the functions "near" and "between" for compiled modules as well as in python.
These are commonly used functions in subdomain definitions, as
near(x[0], coord) -> abs(x[0]-coord) < DOLFIN_EPS,
or
between(coord1, x[0], coord2) -> coord0-DOLFIN_EPS < x[0] < coord1+DOLFIN_EPS.
---
dolfin/math/basic.cpp | 10 +++++++
dolfin/math/basic.h | 5 ++++
site-packages/dolfin/common/__init__.py | 1 +
site-packages/dolfin/common/math.py | 15 +++++++++++
.../dolfin/compilemodules/compilemodules.py | 26 +++++++++++++------
5 files changed, 49 insertions(+), 8 deletions(-)
create mode 100644 site-packages/dolfin/common/math.py
diff --git a/dolfin/math/basic.cpp b/dolfin/math/basic.cpp
index f9215d0..291569e 100644
--- a/dolfin/math/basic.cpp
+++ b/dolfin/math/basic.cpp
@@ -9,6 +9,7 @@
#include <time.h>
#include <cstdlib>
#include <cmath>
+#include <dolfin/common/constants.h>
#include "basic.h"
using namespace dolfin;
@@ -51,3 +52,12 @@ void dolfin::seed(unsigned int s)
rand_seeded = true;
}
//-----------------------------------------------------------------------------
+bool dolfin::near(double x, double x0)
+{
+ return dolfin::between(x0, x, x0);
+}
+//-----------------------------------------------------------------------------
+bool dolfin::between(double x0, double x, double x1)
+{
+ return (x0-DOLFIN_EPS < x && x < x1+DOLFIN_EPS);
+}
diff --git a/dolfin/math/basic.h b/dolfin/math/basic.h
index f4086e9..ac695bf 100644
--- a/dolfin/math/basic.h
+++ b/dolfin/math/basic.h
@@ -25,6 +25,11 @@ namespace dolfin
/// Seed random number generator
void seed(unsigned int s);
+ /// Return true if x is within DOLFIN_EPS of x0
+ bool near(double x, double x0);
+
+ /// Return true if x is between x0 and x1 (inclusive)
+ bool between(double x0, double x, double x1);
}
#endif
diff --git a/site-packages/dolfin/common/__init__.py b/site-packages/dolfin/common/__init__.py
index 72aaf33..4fa6ab3 100644
--- a/site-packages/dolfin/common/__init__.py
+++ b/site-packages/dolfin/common/__init__.py
@@ -4,3 +4,4 @@ from dolfin.common.constants import *
from dolfin.common.plot import *
from dolfin.common.time import *
from dolfin.common.memory import *
+from dolfin.common.math import *
diff --git a/site-packages/dolfin/common/math.py b/site-packages/dolfin/common/math.py
new file mode 100644
index 0000000..0d7b16c
--- /dev/null
+++ b/site-packages/dolfin/common/math.py
@@ -0,0 +1,15 @@
+"""Convenience functions that are useful in subdomains and expressions. Matches
+some of the functions in dolfin/math/basic.h."""
+
+from constants import DOLFIN_EPS
+
+def between(x0, x, x1):
+ """Returns true if x0-epsilon < x < x1+epsilon"""
+ return x0-DOLFIN_EPS < x < x1+DOLFIN_EPS
+
+def near(x, x0):
+ """Returns true if x0-epsilon < x < x0+epsilon"""
+ return x0-DOLFIN_EPS < x < x0+DOLFIN_EPS
+
+def sqr(x):
+ return x*x
diff --git a/site-packages/dolfin/compilemodules/compilemodules.py b/site-packages/dolfin/compilemodules/compilemodules.py
index 525a6d1..d024815 100644
--- a/site-packages/dolfin/compilemodules/compilemodules.py
+++ b/site-packages/dolfin/compilemodules/compilemodules.py
@@ -30,9 +30,18 @@ _math_builtins = [
"exp", "frexp", "ldexp", "log", "log10", "modf",
"pow", "sqrt", "ceil", "fabs", "floor", "fmod"]
-math_header = "\n// cmath functions\n" + \
- "\n".join("using std::%s;"%mf for mf in _math_builtins) + \
- "\nconst double pi = DOLFIN_PI;\n\n"
+_math_dolfin = [
+ # functions from dolfin::math:
+ "sqr", "ipow", "rand", "near", "between"]
+
+math_header = """
+// cmath functions
+
+%s
+
+const double pi = DOLFIN_PI;
+""" % "\n".join("using std::%s;"%mf for mf in _math_builtins)
+
_cpp_keywords = ["auto","const","double","float","int","short","struct","unsigned",
"break","continue","else","for","long","signed","switch","void",
@@ -40,10 +49,10 @@ _cpp_keywords = ["auto","const","double","float","int","short","struct","unsigne
"char","do","extern","if","return","static","union","while",
"asm","dynamic_cast","namespace","reinterpret_cast","try",
"bool","explicit","new","static_cast","typeid","volatile",
- "catch","false","operator","template","typename",
+ "catch","operator","template","typename",
"class","friend","private","this","using",
"const_cast","inline","public","throw","virtual",
- "delete","mutable","protected","true","wchar_t",
+ "delete","mutable","protected","wchar_t",
"or","and","xor","not"]
def expression_to_code_fragments(expr, defaults, arguments):
@@ -71,8 +80,8 @@ def expression_to_code_fragments(expr, defaults, arguments):
# Remove any variables defined in the arguments list
variables.difference_update(arguments)
- # Remove the builtin math functions from the variables (and 'pi', we use DOLFIN_PI)
- variables.difference_update(_math_builtins+["pi"])
+ # Remove the builtin math functions from the variables
+ variables.difference_update(_math_builtins+_math_dolfin+["pi"])
# Remove the numerals from the variables
numerals = [v for v in variables if v[0] in "0123456789"]
@@ -193,7 +202,8 @@ set the environment variable BOOST_DIR.
instant_kwargs['include_dirs'].append(boost_include_dir)
instant_kwargs['system_headers'] = ["cmath", "iostream","complex",
"stdexcept","numpy/arrayobject.h",
- "dolfin.h","boost/shared_ptr.hpp"]
+ "dolfin.h","boost/shared_ptr.hpp",
+ "dolfin/math/basic.h"]
instant_kwargs['swigargs'] =['-c++','-I.']
instant_kwargs['cppargs'] = []
--
1.7.1
From 656d5e8c02c6d703db269c60fee2a463f03a55fb Mon Sep 17 00:00:00 2001
From: Joachim B Haga <jobh@xxxxxxxxx>
Date: Fri, 14 Jan 2011 13:49:36 +0100
Subject: [PATCH 1/9] Let PETScKrylovSolver report reason for divergence as string.
---
dolfin/la/PETScKrylovSolver.cpp | 5 +++--
1 files changed, 3 insertions(+), 2 deletions(-)
diff --git a/dolfin/la/PETScKrylovSolver.cpp b/dolfin/la/PETScKrylovSolver.cpp
index a03ee14..2fdc1bb 100644
--- a/dolfin/la/PETScKrylovSolver.cpp
+++ b/dolfin/la/PETScKrylovSolver.cpp
@@ -201,11 +201,12 @@ dolfin::uint PETScKrylovSolver::solve(PETScVector& x, const PETScVector& b)
// Get solver residual norm
double rnorm = 0.0;
KSPGetResidualNorm(*_ksp, &rnorm);
+ const char *reason_str = KSPConvergedReasons[reason];
bool error_on_nonconvergence = parameters["error_on_nonconvergence"];
if (error_on_nonconvergence)
- error("PETSc Krylov solver did not converge in %i iterations (PETSc reason %i, norm %e).", num_iterations, reason, rnorm);
+ error("PETSc Krylov solver did not converge in %i iterations (PETSc reason %s, norm %e).", num_iterations, reason_str, rnorm);
else
- warning("Krylov solver did not converge in %i iterations (PETSc reason %i, norm %e).", num_iterations, reason, rnorm);
+ warning("Krylov solver did not converge in %i iterations (PETSc reason %s, norm %e).", num_iterations, reason_str, rnorm);
}
// Report results
--
1.7.1
From 4289178bb6fb66dbd06b05f70c98d7654b5e7232 Mon Sep 17 00:00:00 2001
From: Joachim B Haga <jobh@xxxxxxxxx>
Date: Thu, 20 Jan 2011 10:36:36 +0100
Subject: [PATCH 2/9] Changes required to compile with petsc-dev.
---
cmake/modules/FindPETSc.cmake | 4 +++-
dolfin/la/PETScLUSolver.cpp | 12 ++++++++++++
dolfin/la/PETScMatrix.cpp | 18 ++++++++++++++++--
dolfin/la/PETScVector.cpp | 4 ++++
4 files changed, 35 insertions(+), 3 deletions(-)
diff --git a/cmake/modules/FindPETSc.cmake b/cmake/modules/FindPETSc.cmake
index 6875162..e226851 100644
--- a/cmake/modules/FindPETSc.cmake
+++ b/cmake/modules/FindPETSc.cmake
@@ -111,7 +111,9 @@ show :
endmacro()
# Call macro to get the PETSc variables
- petsc_get_variable(PETSC_INCLUDE PETSC_INCLUDE)
+ petsc_get_variable(PETSC_INCLUDE PETSC_INCLUDE) # 3.1
+ petsc_get_variable(PETSC_CC_INCLUDES PETSC_CC_INCLUDES) # dev
+ set(PETSC_INCLUDE ${PETSC_INCLUDE} ${PETSC_CC_INCLUDES})
petsc_get_variable(PETSC_LIB PETSC_LIB)
# Remove temporary Makefile
diff --git a/dolfin/la/PETScLUSolver.cpp b/dolfin/la/PETScLUSolver.cpp
index 142123e..29e2407 100644
--- a/dolfin/la/PETScLUSolver.cpp
+++ b/dolfin/la/PETScLUSolver.cpp
@@ -38,6 +38,18 @@ namespace dolfin
}
};
}
+
+// Forward compatibility with petsc-dev
+#if (PETSC_VERSION_RELEASE == 0)
+#define MAT_SOLVER_UMFPACK MATSOLVERUMFPACK
+#define MAT_SOLVER_MUMPS MATSOLVERMUMPS
+#define MAT_SOLVER_PASTIX MATSOLVERPASTIX
+#define MAT_SOLVER_PETSC MATSOLVERPETSC
+#define MAT_SOLVER_SPOOLES MATSOLVERSPOOLES
+#define MAT_SOLVER_SUPERLU_DIST MATSOLVERSUPERLU_DIST
+#define MAT_SOLVER_SUPERLU MATSOLVERSUPERLU
+#endif
+
//-----------------------------------------------------------------------------
// Available LU solver
const std::map<std::string, const MatSolverPackage> PETScLUSolver::lu_packages
diff --git a/dolfin/la/PETScMatrix.cpp b/dolfin/la/PETScMatrix.cpp
index 343b49b..facc407 100644
--- a/dolfin/la/PETScMatrix.cpp
+++ b/dolfin/la/PETScMatrix.cpp
@@ -288,10 +288,17 @@ void PETScMatrix::zero(uint m, const uint* rows)
assert(A);
IS is = 0;
+ PetscScalar null = 0.0;
+#if (PETSC_VERSION_RELEASE == 0)
+ ISCreateGeneral(PETSC_COMM_SELF, static_cast<int>(m),
+ reinterpret_cast<const int*>(rows),
+ PETSC_COPY_VALUES, &is);
+ MatZeroRowsIS(*A, is, null, NULL, NULL);
+#else
ISCreateGeneral(PETSC_COMM_SELF, static_cast<int>(m),
reinterpret_cast<const int*>(rows), &is);
- PetscScalar null = 0.0;
MatZeroRowsIS(*A, is, null);
+#endif
ISDestroy(is);
}
//-----------------------------------------------------------------------------
@@ -300,10 +307,17 @@ void PETScMatrix::ident(uint m, const uint* rows)
assert(A);
IS is = 0;
+ PetscScalar one = 1.0;
+#if (PETSC_VERSION_RELEASE == 0)
+ ISCreateGeneral(PETSC_COMM_SELF, static_cast<int>(m),
+ reinterpret_cast<const int*>(rows),
+ PETSC_COPY_VALUES, &is);
+ MatZeroRowsIS(*A, is, one, NULL, NULL);
+#else
ISCreateGeneral(PETSC_COMM_SELF, static_cast<int>(m),
reinterpret_cast<const int*>(rows), &is);
- PetscScalar one = 1.0;
MatZeroRowsIS(*A, is, one);
+#endif
ISDestroy(is);
}
//-----------------------------------------------------------------------------
diff --git a/dolfin/la/PETScVector.cpp b/dolfin/la/PETScVector.cpp
index 2933a29..bccb5a1 100644
--- a/dolfin/la/PETScVector.cpp
+++ b/dolfin/la/PETScVector.cpp
@@ -612,7 +612,11 @@ void PETScVector::gather(GenericVector& y, const Array<uint>& indices) const
// Create local index sets
IS from, to;
+#if (PETSC_VERSION_RELEASE == 0)
+ ISCreateGeneral(PETSC_COMM_SELF, n, global_indices, PETSC_COPY_VALUES, &from);
+#else
ISCreateGeneral(PETSC_COMM_SELF, n, global_indices, &from);
+#endif
ISCreateStride(PETSC_COMM_SELF, n, 0 , 1, &to);
// Resize vector if required
--
1.7.1
From 11d341fcbf22e77bf2a6c359415c890c04bacad1 Mon Sep 17 00:00:00 2001
From: Joachim B Haga <jobh@xxxxxxxxx>
Date: Sun, 16 Jan 2011 14:13:07 +0100
Subject: [PATCH 3/9] NoDeleter doesn't require type information.
---
dolfin/common/Array.h | 4 ++--
dolfin/common/NoDeleter.h | 5 ++---
dolfin/la/PETScKrylovMatrix.cpp | 4 ++--
dolfin/la/PETScUserPreconditioner.cpp | 4 ++--
4 files changed, 8 insertions(+), 9 deletions(-)
diff --git a/dolfin/common/Array.h b/dolfin/common/Array.h
index af36738..cc72f5f 100644
--- a/dolfin/common/Array.h
+++ b/dolfin/common/Array.h
@@ -44,7 +44,7 @@ namespace dolfin
Array(uint N, boost::shared_array<T> x) : _size(N), x(x) {}
/// Construct array from a pointer. Array will not take ownership.
- Array(uint N, T* x) : _size(N), x(boost::shared_array<T>(x, NoDeleter<T>())) {}
+ Array(uint N, T* x) : _size(N), x(boost::shared_array<T>(x, NoDeleter())) {}
/// Destructor
~Array() {}
@@ -61,7 +61,7 @@ namespace dolfin
void update(uint N, T* _x)
{
_size = N;
- x.reset(_x, NoDeleter<T>());
+ x.reset(_x, NoDeleter());
}
/// Return informal string representation (pretty-print).
diff --git a/dolfin/common/NoDeleter.h b/dolfin/common/NoDeleter.h
index 098638d..ac520ae 100644
--- a/dolfin/common/NoDeleter.h
+++ b/dolfin/common/NoDeleter.h
@@ -16,11 +16,10 @@ namespace dolfin
/// NoDeleter is a customised deleter intended for use with smart pointers.
- template <class T>
class NoDeleter
{
public:
- void operator() (T* p) {}
+ void operator() (const void *p) {}
};
/// Helper function to construct shared pointer with NoDeleter with cleaner syntax
@@ -28,7 +27,7 @@ namespace dolfin
template<class T>
boost::shared_ptr<T> reference_to_no_delete_pointer(T& r)
{
- return boost::shared_ptr<T>(&r, NoDeleter<T>());
+ return boost::shared_ptr<T>(&r, NoDeleter());
}
}
diff --git a/dolfin/la/PETScKrylovMatrix.cpp b/dolfin/la/PETScKrylovMatrix.cpp
index 7bba30b..aab9d0d 100644
--- a/dolfin/la/PETScKrylovMatrix.cpp
+++ b/dolfin/la/PETScKrylovMatrix.cpp
@@ -28,8 +28,8 @@ namespace dolfin
int usermult(Mat A, Vec x, Vec y)
{
// Wrap x and y in a shared_ptr
- boost::shared_ptr<Vec> _x(&x, NoDeleter<Vec>());
- boost::shared_ptr<Vec> _y(&y, NoDeleter<Vec>());
+ boost::shared_ptr<Vec> _x(&x, NoDeleter());
+ boost::shared_ptr<Vec> _y(&y, NoDeleter());
void* ctx = 0;
MatShellGetContext(A, &ctx);
diff --git a/dolfin/la/PETScUserPreconditioner.cpp b/dolfin/la/PETScUserPreconditioner.cpp
index 82b1aea..450fb28 100644
--- a/dolfin/la/PETScUserPreconditioner.cpp
+++ b/dolfin/la/PETScUserPreconditioner.cpp
@@ -48,8 +48,8 @@ int PETScUserPreconditioner::PCApply(PC pc, Vec x, Vec y)
PETScUserPreconditioner* newpc = (PETScUserPreconditioner*)pc->data;
- boost::shared_ptr<Vec> _x(&x, NoDeleter<Vec>());
- boost::shared_ptr<Vec> _y(&y, NoDeleter<Vec>());
+ boost::shared_ptr<Vec> _x(&x, NoDeleter());
+ boost::shared_ptr<Vec> _y(&y, NoDeleter());
PETScVector dolfinx(_x), dolfiny(_y);
newpc->solve(dolfiny, dolfinx);
--
1.7.1
From a7b5aeb944523a2358c1f9a23203ca5a12683fe3 Mon Sep 17 00:00:00 2001
From: Joachim B Haga <jobh@xxxxxxxxx>
Date: Wed, 19 Jan 2011 22:03:25 +0100
Subject: [PATCH 4/9] Fix epetra crash with assembly of rectangular matrices.
---
dolfin/la/EpetraMatrix.cpp | 40 +++++++++++++++++++++++++++++++---------
dolfin/la/SparsityPattern.cpp | 2 +-
2 files changed, 32 insertions(+), 10 deletions(-)
diff --git a/dolfin/la/EpetraMatrix.cpp b/dolfin/la/EpetraMatrix.cpp
index 3d6a05e..db7d842 100644
--- a/dolfin/la/EpetraMatrix.cpp
+++ b/dolfin/la/EpetraMatrix.cpp
@@ -92,6 +92,14 @@ void EpetraMatrix::init(const GenericSparsityPattern& sparsity_pattern)
Epetra_MpiComm comm = f.get_mpi_comm();
Epetra_Map row_map(sparsity_pattern.size(0), num_local_rows, 0, comm);
+ // For rectangular matrices with more columns than rows, the columns which are
+ // larger than those in row_map are marked as nonlocal (and assembly fails).
+ // The domain_map fixes that problem, at least in the serial case.
+ // FIXME: Needs attention in the parallel case. Maybe range_map is also req'd.
+ const std::pair<uint,uint> colrange = sparsity_pattern.local_range(1);
+ const int num_local_cols = colrange.second - colrange.first;
+ Epetra_Map domain_map(sparsity_pattern.size(1), num_local_cols, 0, comm);
+
// Create Epetra_FECrsGraph
Epetra_FECrsGraph matrix_map(Copy, row_map, reinterpret_cast<int*>(&dnum_nonzeros[0]));
@@ -105,18 +113,32 @@ void EpetraMatrix::init(const GenericSparsityPattern& sparsity_pattern)
matrix_map.InsertGlobalIndices(global_row, row_set->size(), reinterpret_cast<int*>(&_nz_entries[0]));
}
- // Add off-diagonal block indices
- for (row_set = o_pattern.begin(); row_set != o_pattern.end(); ++row_set)
+ for (uint local_row = 0; local_row < d_pattern.size(); local_row++)
{
- const uint global_row = row_set - o_pattern.begin() + n0;
- const std::vector<dolfin::uint>& nz_entries = *row_set;
- std::vector<dolfin::uint>& _nz_entries = const_cast<std::vector<dolfin::uint>& >(nz_entries);
- matrix_map.InsertGlobalIndices(global_row, row_set->size(), reinterpret_cast<int*>(&_nz_entries[0]));
+ const uint global_row = local_row + n0;
+ std::vector<uint> &entries = const_cast<std::vector<uint>&>(d_pattern[local_row]);
+ matrix_map.InsertGlobalIndices(global_row, entries.size(), reinterpret_cast<int*>(&entries[0]));
}
- // Finalise map
- matrix_map.GlobalAssemble();
- matrix_map.OptimizeStorage();
+ // Add off-diagonal block indices (parallel only)
+ for (uint local_row = 0; local_row < o_pattern.size(); local_row++)
+ {
+ const uint global_row = local_row + n0;
+ std::vector<uint> &entries = const_cast<std::vector<uint>&>(o_pattern[local_row]);
+ matrix_map.InsertGlobalIndices(global_row, entries.size(), reinterpret_cast<int*>(&entries[0]));
+ }
+
+ try
+ {
+ // Finalise map. Here, row_map is standing in for RangeMap, which is
+ // probably ok but should be double-checked.
+ matrix_map.GlobalAssemble(domain_map, row_map);
+ matrix_map.OptimizeStorage();
+ }
+ catch (int err)
+ {
+ error("Epetra threw error %d in assemble", err);
+ }
// Create matrix
A.reset( new Epetra_FECrsMatrix(Copy, matrix_map) );
diff --git a/dolfin/la/SparsityPattern.cpp b/dolfin/la/SparsityPattern.cpp
index 9be303a..ab6eb7a 100644
--- a/dolfin/la/SparsityPattern.cpp
+++ b/dolfin/la/SparsityPattern.cpp
@@ -299,7 +299,7 @@ std::vector<std::vector<dolfin::uint> > SparsityPattern::diagonal_pattern(Type t
//-----------------------------------------------------------------------------
std::vector<std::vector<dolfin::uint> > SparsityPattern::off_diagonal_pattern(Type type) const
{
- std::vector<std::vector<uint> > v(diagonal.size());
+ std::vector<std::vector<uint> > v(off_diagonal.size());
for (uint i = 0; i < off_diagonal.size(); ++i)
v[i] = std::vector<uint>(off_diagonal[i].begin(), off_diagonal[i].end());
--
1.7.1
From e29c96829df24bd2bba2f9641ed576b4c5d6876b Mon Sep 17 00:00:00 2001
From: Joachim B Haga <jobh@xxxxxxxxx>
Date: Thu, 20 Jan 2011 01:20:31 +0100
Subject: [PATCH 5/9] Implement DirichletBC::zero_columns for symmetric modification of off-diagonal block matrices.
It optionally takes a non-zero value to put on the diagonal, which if set
can be used for symmetric modification of diagonal blocks, +1 for PD matrices
and -1 for ND matrices.
---
dolfin/fem/DirichletBC.cpp | 76 ++++++++++++++++++++++++++++++++++++++++++++
dolfin/fem/DirichletBC.h | 5 +++
dolfin/la/EpetraMatrix.cpp | 10 +++++-
3 files changed, 90 insertions(+), 1 deletions(-)
diff --git a/dolfin/fem/DirichletBC.cpp b/dolfin/fem/DirichletBC.cpp
index ad6dc16..e09b311 100644
--- a/dolfin/fem/DirichletBC.cpp
+++ b/dolfin/fem/DirichletBC.cpp
@@ -209,6 +209,82 @@ void DirichletBC::zero(GenericMatrix& A) const
A.apply("insert");
}
//-----------------------------------------------------------------------------
+void DirichletBC::zero_columns(GenericMatrix& A, GenericVector& b, double diag_val) const
+{
+ Map bv_map;
+ get_boundary_values(bv_map, _method);
+
+ // Create lookup table of dofs
+ //const uint nrows = A.size(0); // should be equal to b.size()
+ const uint ncols = A.size(1); // should be equal to max possible dof+1
+
+ std::pair<uint,uint> rows = A.local_range(0);
+ //std::pair<uint,uint> cols = A.local_range(1);
+
+ std::vector<char> is_bc_dof(ncols);
+ std::vector<double> bc_dof_val(ncols);
+ for (Map::const_iterator bv = bv_map.begin(); bv != bv_map.end(); ++bv)
+ {
+ is_bc_dof[bv->first] = 1;
+ bc_dof_val[bv->first] = bv->second;
+ }
+
+ // Scan through all columns of all rows, setting to zero if is_bc_dof[column]
+ // At the same time, we collect corrections to the RHS
+
+ std::vector<uint> cols;
+ std::vector<double> vals;
+ std::vector<double> b_vals;
+ std::vector<uint> b_rows;
+
+ for (uint row=rows.first; row<rows.second; row++)
+ {
+ // If diag_val is nonzero, the matrix is a diagonal block (nrows==ncols),
+ // and we can set the whole BC row
+ if (diag_val != 0.0 && is_bc_dof[row])
+ {
+ A.getrow(row, cols, vals);
+ for (uint j=0; j<cols.size(); j++)
+ vals[j] = (cols[j]==row) * diag_val;
+ A.setrow(row, cols, vals);
+ A.apply("insert");
+ b.setitem(row, bc_dof_val[row]*diag_val);
+ }
+ else // Otherwise, we scan the row for BC columns
+ {
+ A.getrow(row, cols, vals);
+ bool row_changed=false;
+ for (uint j=0; j<cols.size(); j++)
+ {
+ const uint col = cols[j];
+
+ // Skip columns that aren't BC, and entries that are zero
+ if (!is_bc_dof[col] || vals[j] == 0.0)
+ continue;
+
+ // We're going to change the row, so make room for it
+ if (!row_changed)
+ {
+ row_changed = true;
+ b_rows.push_back(row);
+ b_vals.push_back(0.0);
+ }
+
+ b_vals.back() -= bc_dof_val[col]*vals[j];
+ vals[j] = 0.0;
+ }
+ if (row_changed)
+ {
+ A.setrow(row, cols, vals);
+ A.apply("insert");
+ }
+ }
+ }
+
+ b.add(&b_vals.front(), b_rows.size(), &b_rows.front());
+ b.apply("add");
+}
+//-----------------------------------------------------------------------------
const std::vector<std::pair<dolfin::uint, dolfin::uint> >& DirichletBC::markers()
{
return facets;
diff --git a/dolfin/fem/DirichletBC.h b/dolfin/fem/DirichletBC.h
index 6cfcda8..0a58ece 100644
--- a/dolfin/fem/DirichletBC.h
+++ b/dolfin/fem/DirichletBC.h
@@ -169,6 +169,11 @@ namespace dolfin
/// non-diagonal matrices in a block matrix.
void zero(GenericMatrix& A) const;
+ /// Make columns associated with boundary conditions zero, and
+ /// update the RHS to reflect the changes. Useful for non-diagonals.
+ /// The diag_val parameter would normally be -1, 0 or 1.
+ void zero_columns(GenericMatrix& A, GenericVector& b, double diag_val=0) const;
+
/// Return boundary markers (facets stored as pairs of cells and local
/// facet numbers)
const std::vector<std::pair<uint, uint> >& markers();
diff --git a/dolfin/la/EpetraMatrix.cpp b/dolfin/la/EpetraMatrix.cpp
index db7d842..82e749d 100644
--- a/dolfin/la/EpetraMatrix.cpp
+++ b/dolfin/la/EpetraMatrix.cpp
@@ -491,7 +491,15 @@ void EpetraMatrix::getrow(uint row, std::vector<uint>& columns,
void EpetraMatrix::setrow(uint row, const std::vector<uint>& columns,
const std::vector<double>& values)
{
- dolfin_not_implemented();
+ static bool print_msg_once=true;
+ if (print_msg_once)
+ {
+ info("EpetraMatrix::setrow is implemented inefficiently");
+ print_msg_once = false;
+ }
+
+ for (uint i=0; i<columns.size(); i++)
+ set(&values[i], 1, &row, 1, &columns[i]);
}
//-----------------------------------------------------------------------------
LinearAlgebraFactory& EpetraMatrix::factory() const
--
1.7.1
From ed21767db782f574cb25183de5a48cf77b04aa6e Mon Sep 17 00:00:00 2001
From: Joachim B Haga <jobh@xxxxxxxxx>
Date: Wed, 26 Jan 2011 14:44:49 +0100
Subject: [PATCH 6/9] Change all checks for instanceof(x,(int,float)) to numpy.isscalar in python.
In practical terms, this means that for example numpy.float64 is now supported.
---
dolfin/swig/la_post.i | 71 ++++++++++++++++----------
dolfin/swig/parameter_post.i | 12 ++--
site-packages/dolfin/function/expression.py | 26 +++------
3 files changed, 59 insertions(+), 50 deletions(-)
diff --git a/dolfin/swig/la_post.i b/dolfin/swig/la_post.i
index 3bc3808..6d7836c 100644
--- a/dolfin/swig/la_post.i
+++ b/dolfin/swig/la_post.i
@@ -180,47 +180,54 @@ PyObject* _get_eigenpair(dolfin::PETScVector& r, dolfin::PETScVector& c, const i
return v
def __contains__(self, value):
- if not isinstance(value, (int, float)):
+ from numpy import isscalar
+ if not isscalar(value):
raise TypeError, "expected scalar"
return _contains(self,value)
def __gt__(self, value):
- if isinstance(value, (int, float)):
+ from numpy import isscalar
+ if not isscalar(value):
return _compare_vector_with_value(self, value, dolfin_gt)
if isinstance(value, GenericVector):
return _compare_vector_with_vector(self, value, dolfin_gt)
return NotImplemented
def __ge__(self,value):
- if isinstance(value, (int, float)):
+ from numpy import isscalar
+ if isscalar(value):
return _compare_vector_with_value(self, value, dolfin_ge)
if isinstance(value, GenericVector):
return _compare_vector_with_vector(self, value, dolfin_ge)
return NotImplemented
def __lt__(self,value):
- if isinstance(value, (int, float)):
+ from numpy import isscalar
+ if isscalar(value):
return _compare_vector_with_value(self, value, dolfin_lt)
if isinstance(value, GenericVector):
return _compare_vector_with_vector(self, value, dolfin_lt)
return NotImplemented
def __le__(self,value):
- if isinstance(value, (int, float)):
+ from numpy import isscalar
+ if isscalar(value):
return _compare_vector_with_value(self, value, dolfin_le)
if isinstance(value, GenericVector):
return _compare_vector_with_vector(self, value, dolfin_le)
return NotImplemented
def __eq__(self,value):
- if isinstance(value, (int, float)):
+ from numpy import isscalar
+ if isscalar(value):
return _compare_vector_with_value(self, value, dolfin_eq)
if isinstance(value, GenericVector):
return _compare_vector_with_vector(self, value, dolfin_eq)
return NotImplemented
def __neq__(self,value):
- if isinstance(value, (int, float)):
+ from numpy import isscalar
+ if isscalar(value):
return _compare_vector_with_value(self, value, dolfin_neq)
if isinstance(value, GenericVector):
return _compare_vector_with_vector(self, value, dolfin_neq)
@@ -238,12 +245,13 @@ PyObject* _get_eigenpair(dolfin::PETScVector& r, dolfin::PETScVector& c, const i
raise ValueError, "cannot delete Vector elements"
def __setslice__(self, i, j, values):
- if i == 0 and (j >= len(self) or j == -1) and isinstance(values, (float, int, GenericVector)):
- if isinstance(values, (float, int)) or len(values) == len(self):
+ if i == 0 and (j >= len(self) or j == -1): # slice == whole
+ if isinstance(values, GenericVector) and len(values) != len(self):
+ raise ValueError, "dimension error"
+ from numpy import isscalar
+ if isinstance(values, GenericVector) or isscalar(values):
self._assign(values)
return
- else:
- raise ValueError, "dimension error"
self.__setitem__(slice(i, j, 1), values)
def __getslice__(self, i, j):
@@ -262,15 +270,15 @@ PyObject* _get_eigenpair(dolfin::PETScVector& r, dolfin::PETScVector& c, const i
raise TypeError, "expected an int, slice, list or numpy array of integers"
def __setitem__(self, indices, values):
- from numpy import ndarray, integer
+ from numpy import ndarray, integer, isscalar
from types import SliceType
if isinstance(indices, (int, integer)):
- if isinstance(values,(float, int, integer)):
+ if numpy.isscalar(values):
return _set_vector_items_value(self, indices, values)
else:
raise TypeError, "provide a scalar to set single item"
elif isinstance(indices, (SliceType, ndarray, list)):
- if isinstance(values, (float, int)):
+ if isscalar(values):
_set_vector_items_value(self, indices, values)
elif isinstance(values, GenericVector):
_set_vector_items_vector(self, indices, values)
@@ -306,7 +314,8 @@ PyObject* _get_eigenpair(dolfin::PETScVector& r, dolfin::PETScVector& c, const i
def __mul__(self,other):
"""x.__mul__(y) <==> x*y"""
- if isinstance(other, (int, float)):
+ from numpy import isscalar
+ if isscalar(other):
ret = self.copy()
ret._scale(other)
return ret
@@ -318,7 +327,8 @@ PyObject* _get_eigenpair(dolfin::PETScVector& r, dolfin::PETScVector& c, const i
def __div__(self,other):
"""x.__div__(y) <==> x/y"""
- if isinstance(other, (int, float)):
+ from numpy import isscalar
+ if isscalar(other):
ret = self.copy()
ret._scale(1.0 / other)
return ret
@@ -334,7 +344,8 @@ PyObject* _get_eigenpair(dolfin::PETScVector& r, dolfin::PETScVector& c, const i
def __rmul__(self, other):
"""x.__rmul__(y) <==> y*x"""
- if isinstance(other,(int,float)):
+ from numpy import isscalar
+ if isscalar(other):
ret = self.copy()
ret._scale(other)
return ret
@@ -360,7 +371,8 @@ PyObject* _get_eigenpair(dolfin::PETScVector& r, dolfin::PETScVector& c, const i
def __imul__(self, other):
"""x.__imul__(y) <==> x*y"""
- if isinstance(other, (float, int)):
+ from numpy import isscalar
+ if isscalar(other):
self._scale(other)
return self
if isinstance(other, GenericVector):
@@ -370,7 +382,8 @@ PyObject* _get_eigenpair(dolfin::PETScVector& r, dolfin::PETScVector& c, const i
def __idiv__(self, other):
"""x.__idiv__(y) <==> x/y"""
- if isinstance(other, (float, int)):
+ from numpy import isscalar
+ if isscalar(other):
self._scale(1.0 / other)
return self
return NotImplemented
@@ -467,7 +480,7 @@ PyObject* _get_eigenpair(dolfin::PETScVector& r, dolfin::PETScVector& c, const i
return down_cast(_get_matrix_sub_matrix(self, indices[0], indices[1]))
def __setitem__(self, indices, values):
- from numpy import ndarray
+ from numpy import ndarray, isscalar
from types import SliceType
if not (isinstance(indices, tuple) and len(indices) == 2):
raise TypeError, "expected two indices"
@@ -476,7 +489,7 @@ PyObject* _get_eigenpair(dolfin::PETScVector& r, dolfin::PETScVector& c, const i
if isinstance(indices[0], int):
if isinstance(indices[1], int):
- if not isinstance(values, (float, int)):
+ if not isscalar(values):
raise TypeError, "expected scalar for single value assigment"
_set_matrix_single_item(self, indices[0], indices[1], values)
else:
@@ -540,8 +553,8 @@ PyObject* _get_eigenpair(dolfin::PETScVector& r, dolfin::PETScVector& c, const i
def __mul__(self,other):
"""x.__mul__(y) <==> x*y"""
- from numpy import ndarray
- if isinstance(other,(int,float)):
+ from numpy import ndarray, isscalar
+ if isscalar(other):
ret = self.copy()
ret._scale(other)
return ret
@@ -574,7 +587,8 @@ PyObject* _get_eigenpair(dolfin::PETScVector& r, dolfin::PETScVector& c, const i
def __div__(self,other):
"""x.__div__(y) <==> x/y"""
- if isinstance(other,(int,float)):
+ from numpy import isscalar
+ if isscalar(other):
ret = self.copy()
ret._scale(1.0/other)
return ret
@@ -590,7 +604,8 @@ PyObject* _get_eigenpair(dolfin::PETScVector& r, dolfin::PETScVector& c, const i
def __rmul__(self,other):
"""x.__rmul__(y) <==> y*x"""
- if isinstance(other, (int, float)):
+ from numpy import isscalar
+ if isscalar(other):
ret = self.copy()
ret._scale(other)
return ret
@@ -616,14 +631,16 @@ PyObject* _get_eigenpair(dolfin::PETScVector& r, dolfin::PETScVector& c, const i
def __imul__(self,other):
"""x.__imul__(y) <==> x*y"""
- if isinstance(other, (float, int)):
+ from numpy import isscalar
+ if isscalar(other):
self._scale(other)
return self
return NotImplemented
def __idiv__(self,other):
"""x.__idiv__(y) <==> x/y"""
- if isinstance(other, (float, int)):
+ from numpy import isscalar
+ if isscalar(other):
self._scale(1.0 / other)
return self
return NotImplemented
diff --git a/dolfin/swig/parameter_post.i b/dolfin/swig/parameter_post.i
index 95e3606..b4e844e 100644
--- a/dolfin/swig/parameter_post.i
+++ b/dolfin/swig/parameter_post.i
@@ -274,17 +274,17 @@ def __new_Parameter_init__(self,*args,**kwargs):
if len(kwargs) == 0:
return
+ from numpy import isscalar
for key, value in kwargs.iteritems():
if isinstance(value,type(self)):
self.add(value)
elif isinstance(value,tuple):
- if len(value) > 0 and ((isinstance(value[0],str) and len(value) == 2) or \
- (isinstance(value[0],(int,float)) and len(value) == 3)):
- if isinstance(value[0],(float,int)) or \
- (isinstance(value[0],str) and isinstance(value[1],list)):
- self.add(key,*value)
- else:
+ if isscalar(value[0]) and len(value) == 3:
+ self.add(key, *value)
+ elif instanceof(value[0], str) and len(value) == 2:
+ if not isinstance(value[1], list):
raise TypeError, "expected a list as second item of tuple, when first is a 'str'"
+ self.add(key, *value)
else:
raise TypeError,"expected a range tuple of size 2 for 'str' values and 3 for scalars"
else:
diff --git a/site-packages/dolfin/function/expression.py b/site-packages/dolfin/function/expression.py
index 8c3a345..500755a 100644
--- a/site-packages/dolfin/function/expression.py
+++ b/site-packages/dolfin/function/expression.py
@@ -54,6 +54,7 @@ import types
# Import UFL and SWIG-generated extension module (DOLFIN C++)
import ufl
import dolfin.cpp as cpp
+import numpy
# Local imports
from dolfin.compilemodules.expressions import compile_expressions
@@ -289,7 +290,6 @@ def expression__call__(self, *args, **kwargs):
... fv(x[i:i+3], values = values[i:i+3])
"""
- import numpy
if len(args)==0:
raise TypeError, "expected at least 1 argument"
@@ -321,24 +321,16 @@ def expression__call__(self, *args, **kwargs):
# Assume all args are x argument
x = args
- # If only one x argument has been provided
+ # If only one x argument has been provided, unpack it if it's an iterable
if len(x) == 1:
- # Check coordinate argument
- if not isinstance(x[0], (int, float, numpy.ndarray, list, tuple)):
- raise TypeError, "expected a scalar or an iterable as coordinate argument"
- # Check for scalar x
- if isinstance(x[0], (int, float)):
- x = numpy.fromiter(x, 'd')
- else:
+ if hasattr(x[0], '__iter__'):
x = x[0]
- if isinstance(x, (list, tuple)):
- x = numpy.fromiter(x, 'd')
- # If several x arguments have been provided
- else:
- if not all(isinstance(v, (int, float)) for v in x):
- raise TypeError, "expected different number of scalar arguments for the coordinates"
- x = numpy.fromiter(x,'d')
+ # Convert it to an 1D numpy array
+ try:
+ x = numpy.fromiter(x, 'd')
+ except (TypeError, ValueError, AssertionError):
+ raise TypeError, "expected scalar arguments for the coordinates"
if len(x) == 0:
raise TypeError, "coordinate argument too short"
@@ -652,7 +644,7 @@ def _check_defaults(defaults):
for key, val in defaults.iteritems():
if not isinstance(key, str):
raise TypeError, "All keys in 'defaults' must be a 'str'."
- if not isinstance(val, (int, float)):
+ if not numpy.isscalar(val):
raise TypeError, "All values in 'defaults' must be scalars."
def _is_complex_expression(cppcode):
--
1.7.1
From 0d27de20e46105f052d4dbb84433606e4732151f Mon Sep 17 00:00:00 2001
From: Joachim B Haga <jobh@xxxxxxxxx>
Date: Thu, 27 Jan 2011 21:20:27 +0100
Subject: [PATCH 8/9] Support for parallel matrix multiplication (PETSc and Epetra).
---
dolfin/la/EpetraMatrix.cpp | 23 ++++++++++++++++-------
dolfin/la/GenericMatrix.h | 8 ++++++--
dolfin/la/PETScBaseMatrix.h | 16 +++++++++-------
dolfin/la/PETScMatrix.cpp | 18 ++++++++++++++++--
dolfin/swig/la_post.i | 12 +++++++-----
5 files changed, 54 insertions(+), 23 deletions(-)
diff --git a/dolfin/la/EpetraMatrix.cpp b/dolfin/la/EpetraMatrix.cpp
index 82e749d..1329134 100644
--- a/dolfin/la/EpetraMatrix.cpp
+++ b/dolfin/la/EpetraMatrix.cpp
@@ -160,13 +160,16 @@ dolfin::uint EpetraMatrix::size(uint dim) const
std::pair<dolfin::uint, dolfin::uint> EpetraMatrix::local_range(uint dim) const
{
assert(dim < 2);
- if (dim == 1)
- error("Cannot compute columns range for Epetra matrices.");
- const Epetra_BlockMap& row_map = A->RowMap();
- assert(row_map.LinearMap());
+ const Epetra_BlockMap *map;
+ if (dim == 0)
+ map = &A->RangeMap(); // Or RowMap, these will usually be equal.
+ else
+ map = &A->DomainMap(); // NOT ColMap, which is not generally linear.
+
+ assert(map->LinearMap());
- return std::make_pair(row_map.MinMyGID(), row_map.MaxMyGID() + 1);
+ return std::make_pair(map->MinMyGID(), map->MaxMyGID() + 1);
}
//-----------------------------------------------------------------------------
void EpetraMatrix::get(double* block, uint m, const uint* rows,
@@ -420,10 +423,15 @@ void EpetraMatrix::mult(const GenericVector& x_, GenericVector& Ax_) const
if (!Ax)
error("EpetraMatrix::mult: The vector Ax should be of type EpetraVector.");
+ // Probably better to check the domainmap/rangemap for compatibility, than to
+ // check the global sizes. However, this check matches the one in
+ // PETScMatrix.
+
if (size(1) != x->size())
error("EpetraMatrix::mult: Matrix and vector dimensions don't match for matrix-vector product.");
- Ax->resize(size(0));
+ if (size(0) != Ax->size())
+ Ax->resize(local_range(0));
const int err = A->Multiply(false, *(x->vec()), *(Ax->vec()));
if (err != 0)
@@ -444,7 +452,8 @@ void EpetraMatrix::transpmult(const GenericVector& x_, GenericVector& Ax_) const
if (size(0) != x->size())
error("EpetraMatrix::transpmult: Matrix and vector dimensions don't match for (transposed) matrix-vector product.");
- Ax->resize(size(1));
+ if (size(1) != Ax->size())
+ Ax->resize(local_range(1));
const int err = A->Multiply(true, *(x->vec()), *(Ax->vec()));
if (err != 0)
diff --git a/dolfin/la/GenericMatrix.h b/dolfin/la/GenericMatrix.h
index 3ea6a74..943071f 100644
--- a/dolfin/la/GenericMatrix.h
+++ b/dolfin/la/GenericMatrix.h
@@ -122,10 +122,14 @@ namespace dolfin
/// Set given rows to identity matrix
virtual void ident(uint m, const uint* rows) = 0;
- /// Matrix-vector product, y = Ax
+ /// Matrix-vector product, y = Ax. The result vector y should not be a copy
+ /// of x, since they may have different parallel distributions. If unsure,
+ /// pass an empty vector which will be resized to fit.
virtual void mult(const GenericVector& x, GenericVector& y) const = 0;
- /// Matrix-vector product, y = A^T x
+ /// Matrix-vector product, y = A^T x. The result vector y should not be a copy
+ /// of x, since they may have different parallel distributions. If unsure,
+ /// pass an empty vector which will be resized to fit.
virtual void transpmult(const GenericVector& x, GenericVector& y) const = 0;
/// Multiply matrix by given number
diff --git a/dolfin/la/PETScBaseMatrix.h b/dolfin/la/PETScBaseMatrix.h
index 829c494..a4f84c5 100644
--- a/dolfin/la/PETScBaseMatrix.h
+++ b/dolfin/la/PETScBaseMatrix.h
@@ -31,6 +31,7 @@ namespace dolfin
}
};
+
/// This class is a base class for matrices that can be used in
/// PETScKrylovSolver.
@@ -68,16 +69,17 @@ namespace dolfin
std::pair<uint, uint> local_range(uint dim) const
{
assert(dim <= 1);
- if (dim == 1)
- error("Cannot compute columns range for PETSc matrices.");
+
+ int m(0), n(0);
if (A)
{
- int m(0), n(0);
- MatGetOwnershipRange(*A, &m, &n);
- return std::make_pair(m, n);
+ if (dim == 0)
+ MatGetOwnershipRange(*A, &m, &n);
+ else
+ MatGetOwnershipRangeColumn(*A, &m, &n);
}
- else
- return std::make_pair(0, 0);
+
+ return std::make_pair(m, n);
}
/// Return PETSc Mat pointer
diff --git a/dolfin/la/PETScMatrix.cpp b/dolfin/la/PETScMatrix.cpp
index facc407..75bff30 100644
--- a/dolfin/la/PETScMatrix.cpp
+++ b/dolfin/la/PETScMatrix.cpp
@@ -328,9 +328,20 @@ void PETScMatrix::mult(const GenericVector& x, GenericVector& y) const
const PETScVector& xx = x.down_cast<PETScVector>();
PETScVector& yy = y.down_cast<PETScVector>();
+ // Check only global sizes for compatibility. This is not enough, since the
+ // distribution may be different for a left-vector and a right-vector (even
+ // if they are the same size for a square matrix). But: A proper
+ // compatibility check requires either extra data structures (like Epetra
+ // blockmaps) or parallel communication, since we must avoid having only some
+ // processors calling resize (this might lead to deadlock if the resize
+ // involves MPI communication).
+
if (size(1) != xx.size())
error("Matrix and vector dimensions don't match for matrix-vector product.");
- yy.resize(size(0));
+
+ if (size(0) != yy.size())
+ yy.resize(local_range(0));
+
MatMult(*A, *xx.vec(), *yy.vec());
}
//-----------------------------------------------------------------------------
@@ -343,7 +354,10 @@ void PETScMatrix::transpmult(const GenericVector& x, GenericVector& y) const
if (size(0) != xx.size())
error("Matrix and vector dimensions don't match for matrix-vector product.");
- yy.resize(size(1));
+
+ if (size(1) != yy.size())
+ yy.resize(local_range(1));
+
MatMultTranspose(*A, *xx.vec(), *yy.vec());
}
//-----------------------------------------------------------------------------
diff --git a/dolfin/swig/la_post.i b/dolfin/swig/la_post.i
index 6d7836c..173f28a 100644
--- a/dolfin/swig/la_post.i
+++ b/dolfin/swig/la_post.i
@@ -564,21 +564,23 @@ PyObject* _get_eigenpair(dolfin::PETScVector& r, dolfin::PETScVector& c, const i
if vector_type not in _matrix_vector_mul_map[matrix_type]:
raise TypeError, "Provide a Vector which can be down_casted to ''"%vector_type.__name__
if type(other) == Vector:
- ret = Vector(self.size(0))
+ ret = Vector()
else:
- ret = vector_type(self.size(0))
+ ret = vector_type()
self.mult(other, ret)
return ret
elif isinstance(other, ndarray):
if len(other.shape) != 1:
raise ValueError, "Provide an 1D NumPy array"
+ first, last = self.local_range(1)
+ local_size = last-first
vec_size = other.shape[0]
if vec_size != self.size(1):
- raise ValueError, "Provide a NumPy array with length %d"%self.size(1)
+ raise ValueError, "Provide a NumPy array with length %d"%local_size
vec_type = _matrix_vector_mul_map[get_tensor_type(self)][0]
- vec = vec_type(vec_size)
+ vec = vec_type(local_size)
vec.set_local(other)
- result_vec = vec.copy()
+ result_vec = vec_type()
self.mult(vec, result_vec)
ret = other.copy()
result_vec.get_local(ret)
--
1.7.1
Follow ups
