dolfin team mailing list archive
-
dolfin team
-
Mailing list archive
-
Message #25526
Fwd: [Branch ~dolfin-core/dolfin/trunk] Rev 6630: Avoid unnessesary resize of result vector for A.mult(x, y)
This change is not safe. For a matrix and a vector to be compatible in
parallel, it is not sufficient to check the global size.
Garth
---------- Forwarded message ----------
From: <noreply@xxxxxxxxxxxxx>
Date: 8 March 2012 20:31
Subject: [Branch ~dolfin-core/dolfin/trunk] Rev 6630: Avoid
unnessesary resize of result vector for A.mult(x, y)
To: Garth Wells <gnw20@xxxxxxxxx>
------------------------------------------------------------
revno: 6630
committer: Johan Hake <hake.dev@xxxxxxxxx>
branch nick: work-trunk
timestamp: Thu 2012-03-08 21:25:41 +0100
message:
Avoid unnessesary resize of result vector for A.mult(x, y)
-- Apply patch from Joachim
modified:
ChangeLog
dolfin/la/EpetraMatrix.cpp
dolfin/la/GenericMatrix.h
dolfin/la/PETScMatrix.cpp
dolfin/la/uBLASMatrix.h
dolfin/swig/la/post.i
test/unit/la/python/KrylovSolver.py
--
lp:dolfin
https://code.launchpad.net/~dolfin-core/dolfin/trunk
Your team DOLFIN Core Team is subscribed to branch lp:dolfin.
To unsubscribe from this branch go to
https://code.launchpad.net/~dolfin-core/dolfin/trunk/+edit-subscription
=== modified file 'ChangeLog'
--- ChangeLog 2012-02-13 21:49:10 +0000
+++ ChangeLog 2012-03-08 20:25:41 +0000
@@ -1,3 +1,4 @@
+ - Avoid unnessesary resize of result vector for A*b
- Major speed-up of mesh topology computation
- Add simple 2D and 3D mesh generation (via CGAL)
- Add creation of mesh from triangulations of points (via CGAL)
=== modified file 'dolfin/la/EpetraMatrix.cpp'
--- dolfin/la/EpetraMatrix.cpp 2012-02-24 12:24:53 +0000
+++ dolfin/la/EpetraMatrix.cpp 2012-03-08 20:25:41 +0000
@@ -568,7 +568,17 @@
}
// Resize RHS
- this->resize(Ax, 0);
+ if (Ax.size() == 0)
+ {
+ this->resize(Ax, 0);
+ }
+
+ if (Ax.size() != size(0))
+ {
+ dolfin_error("EpetraMatrix.cpp",
+ "compute matrix-vector product with Epetra matrix",
+ "Vector for matrix-vector result has wrong size");
+ }
dolfin_assert(x.vec());
dolfin_assert(Ax.vec());
@@ -595,7 +605,17 @@
}
// Resize RHS
- this->resize(Ax, 1);
+ if (Ax.size() == 0)
+ {
+ this->resize(Ax, 1);
+ }
+
+ if (Ax.size() != size(1))
+ {
+ dolfin_error("EpetraMatrix.cpp",
+ "compute transpose matrix-vector product with Epetra matrix",
+ "Vector for transpose matrix-vector result has wrong size");
+ }
const int err = A->Multiply(true, *(x.vec()), *(Ax.vec()));
if (err != 0)
=== modified file 'dolfin/la/GenericMatrix.h'
--- dolfin/la/GenericMatrix.h 2012-02-28 22:23:00 +0000
+++ dolfin/la/GenericMatrix.h 2012-03-08 20:25:41 +0000
@@ -137,10 +137,12 @@
/// 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 y vector must either be zero-sized
+ /// or have correct size and parallel layout.
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 y vector must either be
+ /// zero-sized or have correct size and parallel layout.
virtual void transpmult(const GenericVector& x, GenericVector& y) const = 0;
/// Multiply matrix by given number
=== modified file 'dolfin/la/PETScMatrix.cpp'
--- dolfin/la/PETScMatrix.cpp 2012-02-24 12:24:53 +0000
+++ dolfin/la/PETScMatrix.cpp 2012-03-08 20:25:41 +0000
@@ -393,8 +393,19 @@
"Non-matching dimensions for matrix-vector product");
}
- // Resize if dimensions does not match
- resize(yy, 0);
+ // Resize RHS if empty
+ if (yy.size() == 0)
+ {
+ resize(yy, 0);
+ }
+
+ if (size(0) != yy.size())
+ {
+ dolfin_error("PETScMatrix.cpp",
+ "compute matrix-vector product with PETSc matrix",
+ "Vector for matrix-vector result has wrong size");
+ }
+
MatMult(*A, *xx.vec(), *yy.vec());
}
//-----------------------------------------------------------------------------
@@ -412,7 +423,19 @@
"Non-matching dimensions for transpose matrix-vector product");
}
- resize(yy, 1);
+ // Resize RHS if empty
+ if (yy.size() == 0)
+ {
+ resize(yy, 1);
+ }
+
+ if (size(1) != yy.size())
+ {
+ dolfin_error("PETScMatrix.cpp",
+ "compute transpose matrix-vector product with PETSc matrix",
+ "Vector for transpose matrix-vector result has wrong size");
+ }
+
MatMultTranspose(*A, *xx.vec(), *yy.vec());
}
//-----------------------------------------------------------------------------
=== modified file 'dolfin/la/uBLASMatrix.h'
--- dolfin/la/uBLASMatrix.h 2012-02-24 12:24:53 +0000
+++ dolfin/la/uBLASMatrix.h 2012-03-08 20:25:41 +0000
@@ -461,7 +461,31 @@
template <typename Mat>
void uBLASMatrix<Mat>::mult(const GenericVector& x, GenericVector& y) const
{
- ublas::axpy_prod(A, x.down_cast<uBLASVector>().vec(),
y.down_cast<uBLASVector>().vec(), true);
+
+ const uBLASVector& xx = x.down_cast<uBLASVector>();
+ uBLASVector& yy = y.down_cast<uBLASVector>();
+
+ if (size(1) != xx.size())
+ {
+ dolfin_error("uBLASMatrix.cpp",
+ "compute matrix-vector product with uBLAS matrix",
+ "Non-matching dimensions for matrix-vector product");
+ }
+
+ // Resize RHS if empty
+ if (yy.size() == 0)
+ {
+ resize(yy, 0);
+ }
+
+ if (size(0) != yy.size())
+ {
+ dolfin_error("uBLASMatrix.cpp",
+ "compute matrix-vector product with uBLAS matrix",
+ "Vector for matrix-vector result has wrong size");
+ }
+
+ ublas::axpy_prod(A, xx.vec(), yy.vec(), true);
}
//-----------------------------------------------------------------------------
template <typename Mat>
=== modified file 'dolfin/swig/la/post.i'
--- dolfin/swig/la/post.i 2012-02-21 21:45:42 +0000
+++ dolfin/swig/la/post.i 2012-03-08 20:25:41 +0000
@@ -556,9 +556,9 @@
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):
=== modified file 'test/unit/la/python/KrylovSolver.py'
--- test/unit/la/python/KrylovSolver.py 2012-02-29 18:44:36 +0000
+++ test/unit/la/python/KrylovSolver.py 2012-03-08 20:25:41 +0000
@@ -84,15 +84,12 @@
def __init__(self, A):
self.A = A
- # Use tmp vector as A.mult will resize the vector
- self.y_tmp = PETScVector()
PETScKrylovMatrix.__init__(self, A.size(0), A.size(1))
def mult(self, x, y):
# Make ordinary matrix vector product
- self.A.mult(x, self.y_tmp)
- y[:] = self.y_tmp
+ self.A.mult(x, y)
class MatrixFreeKrylovMatrix(PETScKrylovMatrix) :
=== modified file 'ChangeLog'
--- ChangeLog 2012-02-13 21:49:10 +0000
+++ ChangeLog 2012-03-08 20:25:41 +0000
@@ -1,3 +1,4 @@
+ - Avoid unnessesary resize of result vector for A*b
- Major speed-up of mesh topology computation
- Add simple 2D and 3D mesh generation (via CGAL)
- Add creation of mesh from triangulations of points (via CGAL)
=== modified file 'dolfin/la/EpetraMatrix.cpp'
--- dolfin/la/EpetraMatrix.cpp 2012-02-24 12:24:53 +0000
+++ dolfin/la/EpetraMatrix.cpp 2012-03-08 20:25:41 +0000
@@ -568,7 +568,17 @@
}
// Resize RHS
- this->resize(Ax, 0);
+ if (Ax.size() == 0)
+ {
+ this->resize(Ax, 0);
+ }
+
+ if (Ax.size() != size(0))
+ {
+ dolfin_error("EpetraMatrix.cpp",
+ "compute matrix-vector product with Epetra matrix",
+ "Vector for matrix-vector result has wrong size");
+ }
dolfin_assert(x.vec());
dolfin_assert(Ax.vec());
@@ -595,7 +605,17 @@
}
// Resize RHS
- this->resize(Ax, 1);
+ if (Ax.size() == 0)
+ {
+ this->resize(Ax, 1);
+ }
+
+ if (Ax.size() != size(1))
+ {
+ dolfin_error("EpetraMatrix.cpp",
+ "compute transpose matrix-vector product with Epetra matrix",
+ "Vector for transpose matrix-vector result has wrong size");
+ }
const int err = A->Multiply(true, *(x.vec()), *(Ax.vec()));
if (err != 0)
=== modified file 'dolfin/la/GenericMatrix.h'
--- dolfin/la/GenericMatrix.h 2012-02-28 22:23:00 +0000
+++ dolfin/la/GenericMatrix.h 2012-03-08 20:25:41 +0000
@@ -137,10 +137,12 @@
/// 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 y vector must either be zero-sized
+ /// or have correct size and parallel layout.
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 y vector must either be
+ /// zero-sized or have correct size and parallel layout.
virtual void transpmult(const GenericVector& x, GenericVector& y) const = 0;
/// Multiply matrix by given number
=== modified file 'dolfin/la/PETScMatrix.cpp'
--- dolfin/la/PETScMatrix.cpp 2012-02-24 12:24:53 +0000
+++ dolfin/la/PETScMatrix.cpp 2012-03-08 20:25:41 +0000
@@ -393,8 +393,19 @@
"Non-matching dimensions for matrix-vector product");
}
- // Resize if dimensions does not match
- resize(yy, 0);
+ // Resize RHS if empty
+ if (yy.size() == 0)
+ {
+ resize(yy, 0);
+ }
+
+ if (size(0) != yy.size())
+ {
+ dolfin_error("PETScMatrix.cpp",
+ "compute matrix-vector product with PETSc matrix",
+ "Vector for matrix-vector result has wrong size");
+ }
+
MatMult(*A, *xx.vec(), *yy.vec());
}
//-----------------------------------------------------------------------------
@@ -412,7 +423,19 @@
"Non-matching dimensions for transpose matrix-vector product");
}
- resize(yy, 1);
+ // Resize RHS if empty
+ if (yy.size() == 0)
+ {
+ resize(yy, 1);
+ }
+
+ if (size(1) != yy.size())
+ {
+ dolfin_error("PETScMatrix.cpp",
+ "compute transpose matrix-vector product with PETSc matrix",
+ "Vector for transpose matrix-vector result has wrong size");
+ }
+
MatMultTranspose(*A, *xx.vec(), *yy.vec());
}
//-----------------------------------------------------------------------------
=== modified file 'dolfin/la/uBLASMatrix.h'
--- dolfin/la/uBLASMatrix.h 2012-02-24 12:24:53 +0000
+++ dolfin/la/uBLASMatrix.h 2012-03-08 20:25:41 +0000
@@ -461,7 +461,31 @@
template <typename Mat>
void uBLASMatrix<Mat>::mult(const GenericVector& x, GenericVector& y) const
{
- ublas::axpy_prod(A, x.down_cast<uBLASVector>().vec(), y.down_cast<uBLASVector>().vec(), true);
+
+ const uBLASVector& xx = x.down_cast<uBLASVector>();
+ uBLASVector& yy = y.down_cast<uBLASVector>();
+
+ if (size(1) != xx.size())
+ {
+ dolfin_error("uBLASMatrix.cpp",
+ "compute matrix-vector product with uBLAS matrix",
+ "Non-matching dimensions for matrix-vector product");
+ }
+
+ // Resize RHS if empty
+ if (yy.size() == 0)
+ {
+ resize(yy, 0);
+ }
+
+ if (size(0) != yy.size())
+ {
+ dolfin_error("uBLASMatrix.cpp",
+ "compute matrix-vector product with uBLAS matrix",
+ "Vector for matrix-vector result has wrong size");
+ }
+
+ ublas::axpy_prod(A, xx.vec(), yy.vec(), true);
}
//-----------------------------------------------------------------------------
template <typename Mat>
=== modified file 'dolfin/swig/la/post.i'
--- dolfin/swig/la/post.i 2012-02-21 21:45:42 +0000
+++ dolfin/swig/la/post.i 2012-03-08 20:25:41 +0000
@@ -556,9 +556,9 @@
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):
=== modified file 'test/unit/la/python/KrylovSolver.py'
--- test/unit/la/python/KrylovSolver.py 2012-02-29 18:44:36 +0000
+++ test/unit/la/python/KrylovSolver.py 2012-03-08 20:25:41 +0000
@@ -84,15 +84,12 @@
def __init__(self, A):
self.A = A
- # Use tmp vector as A.mult will resize the vector
- self.y_tmp = PETScVector()
PETScKrylovMatrix.__init__(self, A.size(0), A.size(1))
def mult(self, x, y):
# Make ordinary matrix vector product
- self.A.mult(x, self.y_tmp)
- y[:] = self.y_tmp
+ self.A.mult(x, y)
class MatrixFreeKrylovMatrix(PETScKrylovMatrix) :
