dolfin team mailing list archive
-
dolfin team
-
Mailing list archive
-
Message #21065
Re: Some patches for dolfin (resend plus new)
On Fri, Jan 28, 2011 at 02:36:26PM +0100, Joachim Berdal Haga wrote:
> 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:
Just to check: the attached patches are the ones that have not yet
been applied?
> 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.
Sounds good.
--
Anders
> -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'] = []
> 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
> 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
> 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);
> 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());
>
> 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
> 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):
> 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)
> _______________________________________________
> Mailing list: https://launchpad.net/~dolfin
> Post to : dolfin@xxxxxxxxxxxxxxxxxxx
> Unsubscribe : https://launchpad.net/~dolfin
> More help : https://help.launchpad.net/ListHelp
Follow ups
References
