On 8 Mar 2012, at 21:43, Johan Hake <hake....@gmail.com> wrote:
> On Thu, Mar 8, 2012 at 9:40 PM, Garth N. Wells <gn...@cam.ac.uk> wrote: >> 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. > > True. The logic is that a user need to provide a result vector with > correct size and parallell structure. If not just pass an empty > vector. The problem with always resizing the result is that one loose > performance for iterative solvers using A.mult(x, y). The only proper solution that I see is that a user is responsible for making sure that the size is right and that the function do no resizing. The change is a half-way house, so I think that it should be reverted. The old code was safe, and the proper, efficient alternative is that no resizing be done. Garth > > Johan > >> Garth >> >> >> >> >> ---------- Forwarded message ---------- >> From: <nore...@launchpad.net> >> 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 <gn...@cam.ac.uk> >> >> >> ------------------------------------------------------------ >> revno: 6630 >> committer: Johan Hake <hake....@gmail.com> >> 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) : >> >> _______________________________________________ >> Mailing list: https://launchpad.net/~dolfin >> Post to : dolfin@lists.launchpad.net >> Unsubscribe : https://launchpad.net/~dolfin >> More help : https://help.launchpad.net/ListHelp >> _______________________________________________ Mailing list: https://launchpad.net/~dolfin Post to : dolfin@lists.launchpad.net Unsubscribe : https://launchpad.net/~dolfin More help : https://help.launchpad.net/ListHelp
